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ABSTRACT 

The secular evolution of an infinitely thin tepid isolated galactic disc made of a finite number of particles is described using the 
inhomogeneous Balescu-Lenard equation. Assuming that only tightly wound transient spirals are present in the disc, a WKB approx¬ 
imation provides a simple and tractable quadrature for the corresponding drift and diffusion coefficients. It provides insight into the 
physical processes at work during the secular diffusion of a self-gravitating discrete disc and makes quantitative predictions on the 
initial variations of the distribution function in action space. 

When applied to the secular evolution of an isolated stationary self-gravitating Mestel disc, this formalism predicts initially the impor¬ 
tance of the corotation resonance in the inner regions of the disc leading to a regime involving radial migration and heating. It predicts 
in particular the formation of a “ridge like” feature in action space, in agreement with simulations, but over-estimates the timescale 
involved in its appearance. Swing amplification is likely to resolve this discrepancy. 

In astrophysics, the inhomogeneous Balescu-Lenard equation and its WKB limit may also describe the secular diffusion of giant 
molecular clouds in galactic discs, the secular migration and segregation of planetesimals in proto-planetary discs, or even the long¬ 
term evolution of population of stars within the Galactic center. 
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1. Introduction 


Understanding the dynamical evolution of galactic discs over 
cosmic times is a long-standing endeavour. Self-gravitating discs 
are cold dynamical systems, for which rotation represents an im¬ 
portant reservoir of free energy. Fluctuations of the potential in¬ 
duced by discrete (possibly distant) encounters may be strongly 
amplified, while resonances tend to confine and localise their dis¬ 
sipation: such small stimuli can lead to discs spontaneous evolu¬ 
tion to distinct equilibria. Quantifying the relative importance of 
this intrinsically driven evolution w.r.t. that driven by the environ¬ 
ment is of renewed interest now that their cosmological environ¬ 
ment is firmly established in the context of the ACDM paradigm. 
The effect of intrinsic susceptibility on secular timescales can be 
addressed in the context of kinetic theory, which takes explicitly 
into account such interactions. 

The kinetic theory of stellar systems is an old, yet fundamen¬ 
tal, topic in astrophysics.’ It was initiated by Jeans (1929) and 
Chandrasekhar (1942) in the context of 3D stellar systems such 
as elliptical galaxies and globular clusters. The kinetic theory 
of Coulombian plasmas was developed in parallel by Landau 
(1936) and Vlasov (1938). When encounters are neglected be¬ 
tween the particles (stars or electric charges), one gets a purely 
mean field equation (Jeans 1915; Vlasov 1938) called the col¬ 
lisionless Boltzmann equation, or the Vlasov equation. When 


’ For an historical account of the development of kinetic theories in 
astrophysics, plasma physics, and for systems with long-range interac¬ 
tions, see the introductions of the references Chavanis (2013a,h). 


encounters are taken into account, one gets a kinetic equation 
that includes a collision term. In early works, the collision term 
was obtained by assuming that a particle experiences a succes¬ 
sion of independent two-body encounters with the other parti¬ 
cles. The corresponding kinetic equation can be derived either 
from the Boltzmann equation by considering a limit of weak de¬ 
flections (Landau 1936), or directly from the general form of the 
Fokker-Planck equation by evaluating the diffusion and drift co¬ 
efficients in a binary collisions approximation (Chandrasekhar 
1949; Rosenbluth et al. 1957). In the case of neutral plasmas, 
the system is spatially homogeneous, so the distribution func¬ 
tion depends only on the velocity, hence the name kinetic theory. 
By contrast, stellar systems are spatially inhomogeneous, so the 
distribution function depends on position and velocity. In early 
works on stellar dynamics, spatial inhomogeneity was taken into 
account in the advection term (Vlasov) but the collisional term 
was calculated by making a local approximation, as if the system 
were homogeneous. In 3D, the collisional term displays a log¬ 
arithmic divergence at large scales (Jeans 1929; Landau 1936; 
Chandrasekhar 1942). In the case of plasmas, this divergence is 
due to the neglect of collective effects that are responsible for De¬ 
bye shielding. Landau (1936) phenomenologically introduced a 
cut-off at the Debye length to regularize the divergence. 

Later on, Balescu (1960) and Lenard (1960) developed a rig¬ 
orous kinetic theory of plasmas, taking collective effects into 
account, and obtained a kinetic equation, the so-called Balescu- 
Lenard equation, that does not present any divergence at large 
scales. The Debye shielding is taken naturally into account in 
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their treatment through the dielectric function (that is absent 
from the Landau equation). In the case of stellar systems, the 
divergence at large scales is solved by the spatial inhomogeneity 
of the system and its finite extent. One can phenomenologically 
introduce a cut-off at the Jeans scale (Weinberg 1993), i.e. at the 
system’s size, which would correspond to the analogue of the 
Debye length in plasma physics, but this ad hoc treatment is not 
fully satisfactory. Furthermore, it cannot be applied to cold (cen- 
trifugally supported) stellar discs where spatial inhomogeneity is 
more crucial than in 3Z). 

A more fruitful procedure is to write the kinetic equation 
with angle-action variables that are the appropriate variables to 
describe spatially inhomogeneous multi-periodic systems. When 
collective effects are neglected, one obtains the inhomogeneous 
Landau equation (Chavanis 2007, 2013b). When collective ef¬ 
fects are accounted for, one gets the inhomogeneous Balescu- 
Lenard equation (Heyvaerts 2010; Chavanis 2012a). For self- 
gravitating systems, where the interaction is attractive (instead 
of being repulsive as in Coulombian plasmas), collective effects 
are responsible for an anti-shielding which tends to increase the 
effective mass of the stars, hence reducing the relaxation time. 
The Balescu-Lenard equation is valid at the order 1 jN in an ex¬ 
pansion of the dynamics in terms of this small parameter, where 
A»1 is the number of stars. Therefore, it takes finite-A ef¬ 
fects into account and describes the evolution of the system on a 
timescale of the order Nto, where is the dynamical time. For 
times t^Nto, or for A—> H-oo, it reduces to the Vlasov equation 
which ignores distant encounters between stars. Although the ki¬ 
netic theory was initially developed for 3D stellar systems, the 
final form of the inhomogeneous Balescu-Lenard equation also 
applies to stellar discs such as those considered in this paper. 

Indeed, the Balescu-Lenard non-linear equation accounts for 
self-driven orbital secular diffusion of a gravitating system in¬ 
duced by the intrinsic shot noise due to its discreteness and the 
corresponding long range correlations. Even though this kinetic 
equation was first written down more than fifty years ago, it has 
hardly ever been applied in its prime context, but only in various 
limits where it reduces to simpler kinetic equations, as discussed 
above. 

In this paper, we will focus on solving explicitly such an 
equation describing the self-gravitating response of a tepid thin 
disc to its own stochastic fluctuating potential induced by its fi¬ 
nite number of components. In this cool regime, the self-gravity 
of the disc can be tracked via a local WKB-like response, which 
in turn allows us to simplify the a priori 2D formalism to an ef¬ 
fective (non degenerate) ID formalism. We will compare the pre¬ 
diction of the WKB limit to a numerical experiment presented in 
the literature, and discuss its diagnosis power and possible limi¬ 
tations. 

The paper is organized as follows. Section 2 briefly presents 
the content of the inhomogeneous Balescu-Lenard equation. Sec¬ 
tion 3 focuses on razor thin axisymmetric galactic discs within 
the WKB approximation. Section 4 investigates the formation 
of a narrow resonant ridge in an isolated self-gravitating Mes- 
tel disc. Finally, section 5 wraps up. Appendix A provides a 
short sketch of the derivation of the Balescu-Lenard equation. 
Appendix B considers the inhomogeneous Balescu-Lenard equa¬ 
tion without collective effects. Appendix D compares it to other 
similar kinetic equations, and in particular its Fokker-Planck 
limit. 


2. The inhomogeneous Balescu-Lenard equation 

We consider a system made of A particles. We suppose that 
the gravitational background, associated to the Hamiltonian Hq, 
is stationary and integrable, so that we may always remap the 
physical phase-space coordinates (x, v) to the angle-actions coor¬ 
dinates (6,J) (Goldstein 1950; Born 1960; Binney & Tremaine 
2008). We also introduce the intrisinc frequencies of the system 
SI defined as 


si(j) = e = 


dHo 

dj ■ 


( 1 ) 


Along the unperturbed trajectories, the angles 0 are 27r-periodic 
evolving with the frequency SI, whereas the actions J are con¬ 
served. To describe the long-term evolution of such a system, 
one assumes that there are two decoupled timescales; a short 
dynamical timescale and a secular timescale of collisional evo¬ 
lution. We assume that the system is always in a virialized sta¬ 
ble state (i.e. is a stable stationary solution of the Vlasov equa¬ 
tion), so that the distribution function can be written as a quasi¬ 
stationary distribution F-F(J,t). This is a function of the ac¬ 
tions only that slowly evolves in time due to stellar encounters 
(finite-A effects). From Heyvaerts (2010) and Chavanis (2012a) 
(see also Appendix A for a short sketch of the derivation), the 
secular evolution, induced by collisional finite-A effects, of such 
a quasi-stationary distribution function F{J, t) is given by the 
inhomogeneous Balescu-Lenard equation which reads 


— 

dt 



6Yi(mi-Sli-m2-Sl2) 

l®m,,m,(/l,/2,»Iri^l)P 


3 d 


, ( 2 ) 


where d is the dimension of the physical space, and where we 
used the shortened notation Slj = ^^(7,). The r.h.s of equation (2) 
is the Balescu-Lenard operator which encompasses the secular 
diffusion due to collisional effects, see figure 1 . Because it is the 
divergence of a flux, this writing ensures that the total number 
of stars is exactly conserved during the secular diffusion. The 
Dirac delta 6 D(ini-Sli-m 2 -Sl 2 ) is the sharp resonance condi¬ 
tion. One must note that this condition allows to describe non¬ 
trivial gravitational interactions. Indeed, it can cause non-local 
resonances by coupling different regions of action-space Ji and 
J 2 . Even for local resonances (i.e. J\ -J 2 ), it can allow for non¬ 
trivial coupling of oscillations, as soon as mi and m 2 have non¬ 
zero components. The coefficients l/|iD,„,_m,(7i,72, w)p repre¬ 
sent the dressed susceptibilities of the system, for which col¬ 
lective effects have been taken into account^. To deal with the 
resolution of the non-local Poisson equation, following Kalnajs’ 
matrix method (Kalnajs 1976), one has to introduce a com¬ 
plete biorthonormal basis of potentials and densities and 

p''^\x) such that 

= 4-nGp^P^ , r ix [^‘'P\x)T p^‘>\x) = - 6 ^. (3) 


^ In this paper, we are not interested in the initial complex mechanism 
of violent relaxation (Lynden-Bell 1967), during which the system gets 
virialized, since we intend to describe the long-term evolution of an 
already and continuously virialized system. 

^ In the secular timescale limit, the amplification through the propaga¬ 
tion of waves between resonances is assumed to be instantaneous, see 
Appendix A. 
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Fig. 1: Top: (a) a set of two resonant orhits in the inertial frame ; (b) 
in the rotating frame in which they are resonant - here through ILR- 
COR coupling. Bottom: (c) Fluctuations of the distribution function in 
action-space caused by finite-A^ effects showing overdensities for the 
blue and red orbits. The dashed lines correspond to 3 contour levels of 
the intrinsic frequency a> = m Q. respectively associated with the reso¬ 
nance vector nil (gray lines) and m 2 (black lines). The two sets of orbits 
satisfy the resonant condition mi -Hi =m 2 fl 2 , and therefore lead to a 
secular diffusion of the orbital structure of the disc according to equa¬ 
tion (2). Note that the resonant orbits need not be caught in the same 
resonance (mi t m 2 ), be close in position space nor in action space. 


Thanks to this basis, the susceptibility coefficients are given by 

-7 = y [I-M(cu)]^] [<aL1(/2)]* , (4) 

where M is the response matrix given by 

= (27r)''y r ij . (5) 

m ^ 

In this expression, corresponds to the Fourier transform 

in angles of the basis elements where we used the con¬ 

vention that the Fourier transform of a function X{0, J) is given 
by 

n (6) 

Note that equation (2) is therefore a non-linear function of F, 
both explicitly via the products of F on the r.h.s., describing the 
effects of the binary collisions (as in the Boltzmann or Landau 
equations), but also implicitly via equation (5), which encom¬ 
passes the collective effects and is specific to the Balescu-Lenard 
equation. Note also importantly that - in contrast to its counter¬ 
part in plasma physics - equation (2) does not assume that lo¬ 
cal encounters drive secular evolution; resonant interactions of 
correlated possibly distant dressed orbits are sourcing long term 
orbital distortions. 


2.1. Content of the diffusion equation 

One may also rewrite the Balescu-Lenard equation (2) as an 
anisotropic Fokker-Planck equation introducing the associated 
drift and diffusion coefficients. It then reads 


dt 



dF 

m |A,„,(7i)F(7i) H- D„„(7i)nii ■ — 


, (7) 


where (Ji) and D„,^ (Ji) are respectively the anisotropic drift 
and diffusion coefficients associated to a given resonance mi, i.e. 
to a given Fourier mode mi in angles. They both secularly de¬ 
pend on the distribution function F, but this dependence has not 
been explicitly written out to shorten the notations. The drift co¬ 
efficients are given by 


A,,,{J)^-7t(27tYYj\t 


6 Yi{mi-Sli-m2-Sl2) dF 


while the diffusion coefficients are given by 


Dmi (Ji)^n(2n) 




6 r,{mi-ili-m 2 -^ 2 ) 


F{J2)- 


(9) 


The rewriting from equation (7) allows us to discuss some im¬ 
portant properties of such anisotropic diffusion equation. We in¬ 
troduce the total flux, Ttoi, associated with this diffusion, which 
reads 


Tioi = J] m [aM) F{ J) + DM) m~]- 


dj, 


( 10 ) 


As a consequence, the Balescu-Lenard diffusion equation given 
by the expressions (2) and (7) takes the shortened form 

O TT 

— = div(y^tot) ■ (11) 

dt 

We then define as M(t) the mass contained in a volume W' of the 
action-space at time t, so that we have 

M(t)= fdJF(J,t). (12) 

J'V 

Thanks to the divergence theorem, the variation of mass in 'V due 
to secular diffusion corresponds to the flux of particles through 
the boundary S of this volume so that 


^-£riordS^J^fdS(m-n) 


dF 

A,„(J)FiJ)+DM)m- — 

dj 


, (13) 


where n is the exterior pointing normal vector. One can note 
in equation (13) that the contribution from a given resonance 
m takes the form of a preferential diffusion in the direction 
m. This diffusion is therefore anisotropic because it is maxi¬ 
mum for nccm and equal to 0 for n m-0. To emphasize the 
anisotropy of the diffusion, one may use the formalism of the 
slow and fast actions (Lynden-Bell 1979; Earn & Lynden-Bell 
1996). For simplicity, we consider the 2D case. For a given reso¬ 
nance m-(mi, m 2 ), we consider the change of coordinates 


Jm 

\m\ 


if - 


J m^ 

\m\ 


(14) 


where and are respectively the slow and fast actions as¬ 
sociated to the resonance m. Here m-'- corresponds to the direc¬ 
tion perpendicular to the resonance so that m^ — (m 2 , -mi), and 
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\m\ - yfnTm. Thanks to the chain rule, for any function X{J), one 
has 


dX , , dX 

m -= \m\ -— 
dj dJi. 


jL-cst. 


(15) 


Introducing the natural vector basis elements e^^ = m/\m\ and 
e{„-m-‘-/\m\ associated with this change of coordinates, the dif¬ 
fusion flux associated with a resonance m takes the form 




A„U)F{J) + \m\D^{J)^ 

OJm 


(16) 


Such a rewriting illustrates the fact that as soon as only one reso¬ 
nance m dominates the secular evolution, the diffusion flux will 
be aligned with this resonance. Hence one will observe a narrow 
mono-dimensional diffusion in the preferential 7,^-direction. 
During this diffusion, particles will conserve their fast action 7^, 
which can therefore be seen as an adiabatic invariant, whereas 
their slow action 7,^ gets to change. This strong anistropy in the 
diffusion is an essential property of the Balescu-Lenard equa¬ 
tion (7). 


3. Thin tepid discs and their WKB iimit 

When implementing the inhomogeneous Balescu-Lenard equa¬ 
tion, one encounters two main difficulties. The first one is the ex¬ 
plicit construction of a mapping (x, v) {0, J) since the Balescu- 
Lenard drift and diffusion coefficients must be computed us¬ 
ing angle-actions coordinates. The second difficulty arises from 
the non-locality of Poisson’s equation which requires to intro¬ 
duce potential basis elements as in equation (’). One can 
then compute the response matrix from equation (5), which 
must subsequently be inverted. The following step is to compute 
the drift and diffusion coefficients from equations (8) and (9) 
which requires to explicitly deal with the resonance constraint 
mi fli= 0. In the case of a 2D axisymmetric disc, one 
may implement a WKB approximation (Liouville 1837; Toomre 
1964; Kalnajs 1965; Lin & Shu 1966; Palmer et al. 1989) which 
assumes that the diffusion of the system is made of tightly wound 
spirals. Such an assumption has two main consequences. First 
of all, Poisson’s equation becomes local, resulting in a diagonal 
response matrix. Moreover, it also entails that all the resonances 
becomes exactly local, allowing an explicit calculation of the res¬ 
onant constraint We now detail these two 

elements: epicyclic approximation and WKB assumption. 


3.1. Epicyclic approximation 

For a sufficiently cold disc (i.e. a disc where the radial excursions 
of the stars are small), one can explicitly build up a mapping 
{x,v)h^{6,J) thanks to the epicyclic approximation. We intro¬ 
duce the polar coordinates {R, (p) to describe the infinitely thin 
galactic disc, and introduce their associated momenta (/?«, p^). 
As the disc at equilibrium is axisymmetric, the stationary Hamil¬ 
tonian of the system reads 


Ho{R,(f>, Pr,P^) 


1 

2 



+ lAo(f^), 


(17) 


where i/tq is the stationary background potential within the disc. 
Because i/^o is axisymmetric, it does not depend on (p, so that p^ 


is a conserved quantity. We may then define the first action of 
the system, the angular momentum 7^, as 



'd(p P4,-p^- R^(p. 


(18) 


For a given value of 7^, the equation of evolution of R is then 
given by 


R = 


dlj/eS 
' dR 


(19) 


where i/^efr is an effective potential defined as 




( 20 ) 


The heart of the epicyclic approximation is to assume that small 
radial excursions can be approximated as harmonic librations. 
For a given value of 7^, we implicitly introduce the guiding ra¬ 
dius Rg as 


0 = 


5lAeff 


dR 


dipo 

dR 


Rg Rl 


( 21 ) 


Here Rg(J^) is the radius for which stars with an angular mo¬ 
mentum of 7^ are on exactly circular orbits. It is important to 
note that the mapping between Rg and 7^ is bijective and unam¬ 
biguous (up to the sign of 7^). We may then define the two fre¬ 
quencies of evolution: Q.^{Rg) the azimuthal frequency and K{Rg) 
the epicyclic frequency as follows 


K\Rg) = 


1 dp/Q 
Rg dR 


dR'^ 


R, RV 

d^ipo 


3R^ 


Jl 

+ 34 . 

R, K 


( 22 ) 


As the radial oscillations are supposed to be small, one may 
perform a Taylor expansion at first order of the evolution equa¬ 
tion (19) in the neighborhood of the minimum R-Rg so that 
R satisfies the differential equation R-—K^{R-Rg). Hence one 
can note that in this limit the evolution of the radius of a star is 
the one of a harmonic oscillator centered on Rg. Up to an initial 
phase, one has therefore R{t)-Rg h-A cos(Kt), where A is the am¬ 
plitude of the radial oscillations. The associated radial action 7r 
is then given by 



1 

^R pr = 2 


2 


(23) 


For Jr - 0, the orbit is circular. Within the epicyclic approxima¬ 
tion, the frequencies of motion along the action-torii, introduced 
in equation (1) are given by = An impor¬ 

tant dynamical consequence of this approximation is that these 
two frequencies are only function of 7^ and do not depend on 
Jr, so that the resonance constraint mi = 0 becomes 

simpler. Finally, one can explicitly construct the mapping be¬ 
tween {R,(p,pR,pip) and {6i{,6^,Jr,J^) (Lynden-Bell & Kalnajs 
1972; Palmer 1994; Binney & Tremaine 2008), which takes at 
first order the form 


(/? - Rg+A cos ( 0 r ) , 
U = 0, 2^0 A 


K Rn 


sin(0R). 


(24) 
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Thanks to this mapping and the definitions of the actions from 
equations (18) and (23), the epicyclic approximation allows us to 
build up an explicit mapping between the physical phase-space 
coordinates and the angle-actions ones. 

Finally, throughout our calculation, we will assume that the 
stationary distribution function of the disc is a Schwarzschild 
distribution function (or locally isothermal DF) given by 


F(Rg,Jr) 


n^(Rg)mg) <Rg)Jr 

n K{Rg) crl(Rg) crl(Rg) 


(25) 


where 1 .{Rg) is the surface density of the disc and cr^(Rg), which 
varies within the disc, represents the radial velocity dispersion of 
the stars at a given radius. Increasing values of cr^ correspond to 
hotter discs that are therefore more stable. 


3.2. The WKB basis 

As we are considering a 2D case, the potential basis elements 
introduced in equation (3) must be written as tl/^P\R,(f>) in 
the disc polar coordinates and must be orthonormal to the asso¬ 
ciated surface density lf-P\R,(f)). Using a WKB approximation 
amounts to building up local basis elements thanks to which the 
response matrix will become diagonal. 

3.2.1. Definition of the basis elements 

We introduce the basis elements 

^lk,.kM 0) ^ gi(k,<!>+KR) (z?) , (26) 


where the window function Sr„(/?) is defined as 


SrSR) = 


1 

{ncr'^yi'^ 


exp 


(R-Rq? 

2cr^ 


(27) 


The basis elements are indexed by three numbers: is an az¬ 

imuthal number which parametrizes the angular component of 
the basis elements, Rq is the radius position in the disc around 
which the Gaussian window is centered, and kr is the radial 
frequency of the basis element. We also introduced an additional 
parameter cr of scale-separation, which will ensure the biorthog¬ 
onality of the basis elements, as detailed later on. Finally, is 
an amplitude which will be tuned in order to normalize correctly 
the basis elements. Thanks to a somewhat unsual normalization 
of we will ensure that ^ is independent of cr. Figure 2 illus¬ 
trates the radial dependence of the basis elements. Figure 3 illus¬ 
trates the shape of these basis elements in the polar {R, (^)-plane. 
The next steps will be to ensure that these WKB basis elements 
have all the properties required to allow for the computation of 
the dressed susceptibility coefficients introduced in equation (4). 
Therefore, we will successively compute the associated surface 
density elements Y^-h’kr.Ra]^ ensure the biorthogonality of the ba¬ 
sis elements and their correct normalization, and finally compute 
the Fourier transform in angles of the basis elements. 


3.2.2. Associated surface densities 

In order to ensure the biorthogonality of the basis, we will first 
build up the surface densities associated to the potential elements 
introduced in equation (26). We extend the WKB potential in the 
z-direction using the Ansatz 

0, z) = Z(z). (28) 


‘k 



Fig. 2: Two WKB basis elements. Each Gaussian is centered around a 
radius Rq. The typical extension of the Gaussian is given by the decou¬ 
pling scale cr, and they are modulated at the radial frequency U. 



Fig. 3: Two WKB basis elements in the polar (R, 0)-plane. Each basis 
element is located around a central radius Rg, on a region of size cr. 
The winding of the spirals is governed by the radial frequency U. The 
number of azimuthal patterns is given by the index k^, e.g. k^ = l for 
the interior dark gray element, whereas k^ = 2 for the exterior light gray 
one. 


Poisson’s equation in vacuum =0 leads to 

1 


7" 

— -K 

z 


i 1 R-Ro 

--1-2;-T--I- 

krR cr^ kr 


1 


{akrf (krR)^ 


R (a-kr)^ 
R-Ro 1 ^2 


0-2 krJ j 


(29) 


We now explicitly introduce our WKB assumptions. We assume 
that the spirals are tightly wound so that 


krR » I. 


(30) 


Introducing the typical size of the system Rsys, we also addition¬ 
ally suppose that we have 


krCT » 


Rsys 

cr 


(31) 
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Assuming that is of the order of unity, equation (29) becomes 


Z" 

Y 




(32) 


Hence within the WKB limit, the extended potential from equa¬ 
tion (28) takes the form 

0, z) = (j)) , (33) 


where we ensured that for z—> +oo the potential tends to 0, there¬ 
fore introducing a discontinuity for dij/ldz in z = 0. Thanks to 
Gauss theorem, the associated surface density satisfies 


Y.{R,cp) 


1 

4^ 


lim ^ - lim ^ 

z-^o+ dz z-^o- dz 


(34) 


therefore consider a spread cr, central radii Rq, and radial fre¬ 
quencies kr such that 

A/?o » cr » tV ■ (41) 

Akr 

With these constraints, one must necessarily have = k^, kf = kf 
and Rg - R‘q in order to have a non negligible term in equa¬ 
tion (36). It then only remains to explicitly estimate the ampli¬ 
tude of the basis elements. Equation (36) gives 




\kr 


\/n(T- 


if 


dRR exp 


(R-Rof 


= 1. 


(42) 


Thanks to the WKB assumption from equation (41), the integra¬ 
tion can be straightforwardly computed and leads to 


so that we have 


2nG 


(35) 


3.2.3. Biorthogonality condition and normalization 

One must now ensure that the basis elements introduced in equa¬ 
tions (26) and (35) form a biorthogonal basis as assumed in equa¬ 
tion (3). Indeed, it has to satisfy the property 





The r.h.s of this expression takes the form 


(36) 




'ncr^J 


</■ 


X dRRe 


no" 

Hk'’,-kf)R 


(R-Rff 


exp 


2cr^ 


exp 


(R-Rl? 


2 <t^ 


(37) 


The integration on (ft is straightforward and is equal to 2n6 I. In 

order to perform the integration on R, we have to introduce ad¬ 
ditional assumptions to ensure the biorthogonality of the basis. 
The peaks of the Gaussians in equation (37) can be considered 
as separated if ARq=Rq-Rq satisfies the condition 


ARo » cr if R^q^RI- 


(38) 


Under this assumption , the term from equation (37) can be as¬ 
sumed to be non-zero only for R^ -Rq. The r.h.s of equation (36) 
then takes the form 


kl Jit l^r I 

^0 (j 




dRRe 


Hkf-k1)R 


exp 


{R-Rf,f 


(39) 


The integration on R takes the form of a radial Fourier transform 
of a Gaussian of spread cr at the frequency Ak^^kf-k^. It is 
therefore proportional to exp[-(Akr)^/(4/cr)^]. Hence we will 
suppose that the frequency spread Akr satisfies 


Akr » — if k^’r + k1 . 
cr 


(40) 


Under this assumption, the term from equation (39) is non-zero 
only for k^’r-k^. In order to have a biorthogonal basis, one must 

This is an assumption that one might want to lift partially to account 
for weakly non local effects. 



(43) 


3.2.4. Fourier transform in angles 


In order to estimate the susceptibility coefficients and the re¬ 
sponse matrix from equations (4) and (5), one has to be able 
to calculate for the WKB basis elements. Thanks to the 

explicit mapping from equation (24), we have to compute 

^ ^i[krA cosiB.yk,-^ A COS(0fi)) . (44) 

The integration on 6 ^ is straighforward and equal to 2n^^^. Re¬ 
garding the dependence on Or in the complex exponential, we 
may write 

krAcos(0R)-k^ —^■^sin(6»R) ^ Hk^(kr)sm(0R + ePj^), (45) 

K Kg 

where the amplitude Hk^ikr) and the phase shift 0^ are given by 


i/K(k,)-A|k,| Jl- 


2 k(i) 


K krRg 


; = tan * 


K krRg 


2 k(h 


(46) 


For typical galaxies, we have II2<Q.^Ik< 1 (Binney & Tremaine 
2008). Assuming that the azimuthal number k^ is of the order 
of unity, one can use the WKB hypothesis introduced in equa¬ 
tion (30), so that equations (46) can be approximated as 


Hk^ikr) - A \kr\ - \kr\ ; ■ (47) 

Because we have assumed that the disc is tepid, the radial oscilla¬ 
tions are small so that A<&Rg. We may then get rid of the depen- 
dances on A in Sn^iRg+A cos(0r)) by replacing it with Sn^iRg), 
so that the only remaining dependence on A will be in the com¬ 
plex exponentials. To be able to explicitly perform the remaining 
integration on Or in equation (44), we introduce the Bessel func¬ 
tions J'l of the first kind which satisfy the relation 

(48) 

feZ 

We then finally obtain the expression of the Fourier transform in 
angles of the WKB basis elements which reads 

^ [Hm^ikr)] Sr„(R,) . (49) 
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3.3. Estimation of the response matrix 

Thanks to the WKB basis introduced in equation (26), one can 
now explicitly compute the response matrix from equation (5). 
Indeed, we use the expression (49) of the Fourier transform of 
the basis elements, and after simplification of the phase-shift 
terms thanks to the approximation from equation (47), one has 
to evaluate 


m dFIdJ ^k’; ki j(k‘i-k';)R, 

co-m-n 








(50) 


The first step of the calculation is to show that in the WKB limit, 
the response matrix becomes diagonal. One should note that 
the previous expression (50) is similar to equation (37), where 
we discussed the biorthogonality of the WKB basis. In equa¬ 
tion (50), the azimuthal Kronecker symbols ensures that =k^. 
Moreover, because of our WKB assumptions from equation (41) 
on the step distances of the basis elements, the product of the 
two radial Gaussians in Rg imposes that R^ -R^ in order to have 
a non-zero contribution. In order to shorten the notations, we 
temporarily introduce the function h{Rg) defined as 


HRg) 


d70 

d^ 


mdF Id J 

cj-mSl 




ILl-P 


Jm 



which encompasses all the additional radial dependences appear¬ 
ing in equation (50). Thanks to the change of variables Rg, 
the integral on which has to be evaluated in equation (50), 
when estimated for R^ =Rq, is qualitatively of the form 


JdRgh(Rg)e‘^^(^'-'^"^ 


exp 




-Kf 


(51) 


This expression corresponds to a radial Fourier transform IF at 
the frequency Ak^. It can be rewritten as a convolution of two 
radial Fourier transforms so that it becomes 


(51)cx Jdk'nk](k') exp 


(Akr-k'f 

4/cr2 


(52) 


where Akr-k^ -k^. We now rely on the WKB assumption from 
equation (40). If one has Akr + 0, because of the Gaussian from 
equation (52), the contribution from Tih] will come from the re¬ 
gion k' ~ Akr » 1 /cr. We assume that the function h is such that 
its Fourier Transform is limited to the frequency region |k'| < 1 jcr. 
This is consistent with assuming that the properties of the disc 
are radially slowly varying, and this implies that non-zero con¬ 
tributions to the response matrix can only be obtained when 
Akr-kr-kf-O. Therefore, we have shown that within our WKB 
formalism, the response matrix from equation (5) is diagonal. 

In order to shorten the notations, we will denote the matrix 
eigenvalues as 


d[k,p,kr,Rt,]((^) - ■ (53) 

For these diagonal coefficients, the last step is to explicitly com¬ 
pute the integrals over and Jr in equation (50) to obtain the 
expression of the response matrix eigenvalues. We now detail 
this calculation. Thanks to our scale-decoupling approach, we 
may replace the radial Gaussian from equation (51) by a Dirac 


delta duiRg-R^) while paying a careful attention to the correct 
normalization of the Gaussian. Hence we have to evaluate 


d[k.p,kr,Ro]i^) 
(27r)W 


dRr 


^0 m 


mdFIdJ 

CJ-m^il^-nirK 


J 

m 


(54) 


Because of the presence of the azimuthal Kronecker symbol, we 
may drop the sum on m^. The intrinsic frequencies from equa¬ 
tions ( 22 ) allow us to compute 




dRr. 


Ro 


2n^ 


(55) 


Moreover, we assume that the galactic disc is tepid so that 
\dF/dJ^\<i:\dF/dJr\. We may then only keep the term corre¬ 
sponding to a gradient with respect to the radial action F- 
Thanks to the expression of the Schwarzschild distribution func¬ 
tion from equation (25) and the expression of the basis ampitude 
from equation (43), equation (54) becomes after some simple 
algebra 


d-ik.p.kr.Roji^jr) 


27rGS|kr 


k?cr4 


fdJ 

-mrexp[-xJr/o-^] 2 

f^k 

1 

1 r\ tUr 

V K 


(56) 


We may now use the following integration formula (see formula 
(6.615) from Gradshteyn & Ryzhik (2007)) 

X +“ r/i2 


2a 


(57) 


where a > 0,13 > 0,mr e Z and are the modified Bessel func¬ 
tions of the first kind. We apply this formula with a-Klcr^ and 
P - yflkffK. We also introduce the notation 

-HI 


o-t 


so that equation (56) becomes 
27rG2|kr| K 


dlk.p,kr,Ro]i(^) 


X 




-nire ^ Jm,[y] 

U)-ksQ.s-mrK 


(58) 


(59) 


We now define the dimensionless shifted frequency i as 
w-k^Q^ 

i = -. (60) 

K 

Because we have I-m, M-Im, Ix]^ we may rewrite equation (59) 
using the reduction factor (Kalnajs 1965; Lin & Shu 1966) de¬ 
fined as 


-X 

r(s,x)^2{l-s^) — y 


ImM 

l-[i/OTf]^ 


(61) 


As a conclusion, we obtain that within our WKB formalism the 
response matrix M becomes diagonal and in the limit of tepid 
discs reads 


M, 


{k:,K,Rii[ki,kiRi\ 



27rGS|kr| 


T{s,x), 


(62) 


This eigenvalue recovered using the WKB basis introduced 
in equation (26) is in full agreement with the seminal results 
from Kalnajs (1965) and Lin & Shu (1966). In order to handle 
the singularity of the eigenvalue appearing for s-neZ, one adds 
a small imaginary part to the frequency of evaluation, so that 
s-n+irf. Indeed, as long as 77 is small compared to the imaginary 
part of the least damped mode of the disc, adding this complex 
part makes a negligible contribution to the expression of Re(/1). 
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3.4. Estimation of the susceptibiiity coefficients 


where the resonance condition = 0 is given by 


One can now estimate the dressed susceptibility coefficients 
from equation (4). In order to shorten the notations, we will write 
the WKB basis elements introduced in equation (26) as 




(63) 


We have shown previously in equation (62) that within our WKB 
basis, the response matrix M is diagonal. Its eigenvalues will be 
noted as Ap so that we have = S^pAp. Hence the expression (4) 
of the susceptibility coefficients takes the form 


1 

,m2 




1 

l-Apioj) 




Using the expression of the Fourier transformed basis elements 
obtained in equation (49), we obtain 


Z k' I 

6\5 

m, I 


k'l X G 1 


fD,„„^,{JuJ2,u) mtkPR^^l-Ap 


^ m 


HkP 


Sfm'j 


: exp 


/TTCr^ 




(Ri-R^of 


2cr2 


exp 




(R 2 -R^of 


2 (X 


(64) 


where we used the shortened notations Ki-K(Ji), and Ri-Rg(Ji). 
We also used the approximation introduced in equation (47) for 
the values at which the Bessel functions have to be evaluated. 
Thanks to the Kronecker symbols in and we necessarily 
have 

k ‘^, (65) 

so that the sum on can be dropped. 


/(«5) = ^^lI■£^(/^I)-^n2■^^(/^2)■ 


(68) 


The radii R^ therefore correspond to the resonant radii for which 
the resonance condition is satisfied. When writing equation (67), 
we have assumed that the zeros of the resonance function are 
simple, which corresponds to the assumption that for any reso¬ 
nant radius R^, we have ffR'i^ + O. Assuming that 4^ 0, this 
condition can be rewritten as 


'«2 

«2 ^2 


(69) 


Resonance poles are therefore simple as long as the rates of 
change of the two intrinsic frequencies are not in a rational ratio. 
One must note that the Keplerian case for which k-D.^ and the 
harmonic case for which k - 20.^ are in this sense degenerate. It 
can lead to resonant poles of higher multiplicity and would there¬ 
fore require a more involved evaluation of the Balescu-Lenard 
collision operator. In what follows we assume that the potential 
is not degenerate. 

Let us now use the properties of the WKB basis to restrict the 
range of resonant radii, Rj. The expression (64) of the suscepti¬ 
bility coefficients, thanks to the two Gaussians, imposes that the 
relevant resonant radius Rj must necessarily be close to Ri. As 
noted in equation (65), in order to have a non-zero susceptibil¬ 
ity, one also has to satisfy the constraint ni^ - m^. The resonant 
condition which has to be satisfied is therefore given by 


m^Q.^(Ri)+m[K{Ri) - fn^n^(R2)+ni2K(R2) , (70) 

where the distance AR = Ri-R 2 is such that |AR|<(few)cr. Be¬ 
cause the scale-decoupling parameter cr is supposed to be small 
compared to the size of the system, we may approximate the pre¬ 
vious resonant condition as 


3.5. Restriction on the ioci of resonance 


Before proceeding with the evaluation of the susceptibility coef¬ 
ficients obtained in equation (64), let us first emphasize a crucial 
consequence of the WKB basis from equation (26), which is the 
restriction to only exactly local resonances. One can note that 
the expressions (8) and (9) of the drift and diffusion coefficients 
all involve an integration over the mute variable J 2 . For a given 
value of Ji, mi and m 2 , this integration should be seen as a scan 
of the entire action-space, searching for resonant region where 
the constraint mi ili-m2-il2-0 is satisfied. We first recall the 
rule for the composition of a Dirac delta and a function which 
reads 


yeZf 


dpjx-y) 

\f'iy)\ ’ 


( 66 ) 


where Zf-{y\ f(y) - 0), and we have supposed that all the poles 
of / are simple. As noted in equation (22), within the epicyclic 
approximation, the intrinsic frequencies Q = (Q^, k) only depend 
on Rg-RgiJ^) and are independent of Jr. Hence, the resonance 
condition /wi ■ £2i -/M 2 ■ ^2 = 0 only depends on 7^ and is indepen¬ 
dent of J^. Hence if we consider fixed J\, m\ and m 2 , the reso¬ 
nant Dirac delta which has to be studied takes the form 


6 y 2 {R 2 -R’ 2 ) 


Soimi-ni-mrn2)^ Z |Ar;„,.oi| 

RJ|/(t?O=0 iil| 


IK 


(67) 


- dR 


-+m.. 


. Ok 

-m\ 


AR = 


k(Ri). 


(71) 


On the l.h.s of equation (7 1 ), the term within bracket is non-zero, 
because we assumed in equation (69) that the resonant poles are 
simple. Moreover, AR is small, because of our scale-decoupling 
approach. The r.h.s of equation (71) is discrete: it is either zero 
or at least of the order of k{Ri). Because the l.h.s is necessarily 
small, we must have 


R2 — Ri, '^2 “ '”1 • 


(72) 


This result is a crucial consequence of our WKB tightly wound 
spiral assumption, because it implies that only local resonances 
are allowed. In particular this implies that the WKB limit does 
not allow for distant orbits to resonate (through e.g. propagation 
of swing amplified wave packets, see below). Then the sum YjR' 
from equation (67) can be limited to the evaluation in 
Hence within this WKB limit, the susceptibility coefficients from 
equation (64) have to be evaluated only for m 2 - mi and R 2 =Ri, 
so that we have to deal with the expression 


1 

fDmt,mi(Rl,Jl-,Rl,Jr^ X) 



G 1 
k’rRl l-Ap 



Jm 



1 


exp 


(Ri-R^of 


cr^ 


(73) 
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3.6. Asymptotic continuous iimit 

One can note that in equation (73) the susceptibility coefficients 
are still expressed as a discrete sum on the basis index and 
Rq. Our next step is to replace these sums by continuous in¬ 
tegrals. The discrete basis elements are separated by the step 
distances ARq and Akr, which must satisfy the WKB hypothe¬ 
sis detailed in equation (41). We use the Riemann sum formula 
Yjfix)Ax^ Jdxfix), where Ajr controls the distance between the 
basis elements. This transformation is a subtle stage of the calcu¬ 
lation, because one has to consider step distances ARq and Akr, 
which have to simultaneously be large to comply with the WKB 
assumption from equation (41) and small to allow the use of the 
Riemann sum formula. As we are going to transform both the 
sums on k'l and R^, the exact value of the susceptibility coeffi¬ 
cients will depend on our choice for ARq One has to con¬ 
sider the case 


ARq Akr - 27r. 


(74) 


This sampling corresponds to a critical sampling condition 
(Gabor 1946; Daubechies 1990) (See also Fouvry et al. 2015). 
Equation (73) then takes the form 


1 

,mi (R\,JI,R\,J},CU) 


G_ 

2 n 




Akr ARi 


0 


1 1 
krRo l-AkXRo,(^) 


^ AT 





(Ri-Rof 


(75) 


One can now assume that the radial Gaussian present in equa¬ 
tion (75) is sufficiently peaked. Because it is correctly normal¬ 
ized, we may in this limit replace it by (5d(Ri-Ro). The integra¬ 
tion on Rq can then be immediately performed to give 


1 


1 


X 




■A/i (Ri,a>) 


Ztt Ri 

Jm 


Jm 



(76) 


where only A^^ depends on the frequency of evaluation co. One 
may note that in equation (76), all the dependencies in cr have 
disappeared, so that the value of the susceptibility coefficients is 
independent of the precise choice of the WKB basis.The square 
of the susceptibility coefficients which is required to estimate 
the drift and diffusion coefficients from equation (8) and (9) is 
therefore given by 




J_G^ 

4 -n^ r2 


X 


1 1 



(77) 


where we introduced a cut-off at 1/cr^ for the integration on 
kr. This bound is justified by the WKB constraint from equa¬ 
tion (41), which imposes that the probed radial frequency region 
is bounded from below. It is also important to note that these sus¬ 
ceptibility coefficients should be evaluated at since we 

proved in equation (72) that, consistently with our WKB approx¬ 
imation, exactly local resonances are the only ones which have 
to be considered. 

At this stage, there is an arbitration to make between two 
possible behaviors depending on the physical properties of 
the underlying disc. First of all, if the amplification function 
kr>-^l/{l-Ai;A) is asymptotically a sharp function reaching a 


maximum value Tmax for kr-kra^x, one can assume that the sus¬ 
ceptibility coefficients are dominated by the contribution from 
the peak in Au^. In this situation, we can perform an approxima¬ 
tion of the small denominators. The second possible behavior 
arises if the function kri-> 1/(1-/!*•_.) is asymptotically flat, so 
that there is no characteristic scale of blow-up of the amplifi¬ 
cation eigenvalues. In such a situation, the susceptibility coeffi¬ 
cients are mostly dominated by the behavior at the boundaries 
of integration from equation (77) where kr llcri,. The detailed 
response structure of the self-gravitating disc then does not play 
a significant role. 

We place ourselves within the approximation of the small de¬ 
nominators, assuming that the biggest contribution to the suscep¬ 
tibility coefficients comes from waves which yield the largest A),^. 
Therefore, one has to suppose that the function kr^ lUl—AuX 
is a sharp function reaching a maximum value /ln,ax(Ri, w), for 
kr-kmax{R\,<T), with a characteristic spread w). The ex¬ 

pression (77) then becomes 




1 G^(Ak,f 

471-2 p2 k 2 
^71 A j /Cjnax 



(78) 


While still focusing on the contribution to the susceptibility coef¬ 
ficients due to the waves with the largest amplification A{kr), one 
can improve the approximation of the small denominators from 
equation (78). Indeed, starting from equation (77), one can in¬ 
stead perform the integration for k^e [kinf; kjup], where these 
bounds are given by /I(kinf/sup) = dmax/2. This approach is numer¬ 
ically more demanding but does not alter the principal conclu¬ 
sions drawn in this paper, while allowing a more precise determi¬ 
nation of the secular flux structure. All the calculations presented 
in section 4 were performed with this improved approximation. 
Finally, in Appendix B, we detail how this same WKB formalism 
may be applied to the inhomogeneous Balescu-Fenard equation 
without collective effects (Chavanis 2013b). 


3.7. Estimation of the drift and diffusion coefficients 


The drift and diffusion coefficients are given by equations (8) 
and (9). Within the WKB approximation, we have shown in equa¬ 
tions (65) and (72) that the susceptibility coefficients have to be 
evaluated only for nii = m 2 , so that the sum on m 2 in the expres¬ 
sions of the drift and diffusion coefficients may be dropped. As 
the resonances are exactly local, using the formula from equa¬ 
tion (67) adds a prefactor of the form l/\d{mi-£l)/dJ^\, so that 
the drift coefficients from equation (8) become 


Anil (7 1 ) — 


47r^ 


dJ, 


[mi-fii] 




X dJr 


mrdFldJ{j\,J^) 




Similarly, the diffusion coefficients are given by 
47r3 


DmXJl) -7 


\Mmvnx\ 


</dy; 




\D„i,,„iXJl,Jl,Jl,JhmvElX\'^' 


(79) 


(80) 
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In both equations (79) and (80), the susceptibility coefficients are 
given by equation (77) or, within the approximation of the small 
denominators, by equation (78). 

Such simple expressions of the drift and diffusion coeffi¬ 
cients along with the required expression of the susceptibility 
coefficients constitute a main result of this paper. Note impor¬ 
tantly that this WKB formalism is self contained and is ob¬ 
tained without any ad hoc assumptions or fittings. Except for 
the explicit recovery of the amplification eigenvalues from equa¬ 
tion (62), the calculations presented previously are not limited 
to the Schwarzchild distribution function from equation (25). 
Indeed, the drift and diffusion coefficients from equations (79) 
and (80) are valid for any tepid disc, as long as the epicyclic 
angles-actions mapping from equation (24) may be used. In 
Appendix C, we show how the drift and diffusion coefficients 
from equations (79) and (80) can be explicitly computed, when 
one considers a Schwarzschild distribution function as in equa¬ 
tion (25) and when the susceptibility coefficients are estimated 
thanks to the approximation of the small denominators from 
equation (78). Finally, in Appendix D, we compare this 2D WKB 
Balescu-Lenard equation with other similar kinetic equations. 


4. Application 

Let us now illustrate how this WKB approximation of the in¬ 
homogeneous Balescu-Lenard equation can be implemented to 
recover some results obtained via well-crafted numerical simu¬ 
lations. Indeed, Sellwood (2012) (hereafter S12) using careful 
numerical simulations studied the long-term evolution of an iso¬ 
lated stable and stationary Mestel disc (Mestel 1963) sampled 
by pointwise particles. When evolved for hundreds of dynam¬ 
ical times, such a disc would secularly diffuse in action-space 
through the spontaneous generation of transient spiral structures. 
Figure 7 of S12 shows for instance the late time formation of 
resonant ridges along very specific resonant directions. Such fea¬ 
tures are possible signatures of secular evolution, which results 
in a long-term aperiodic evolution of a self-gravitating system, 
during which small resonant and cumulative effects can add up 
in a coherent way. These small effects, which are then amplified 
through the self-gravity of the system originate from finite-N ef¬ 
fects. Indeed, the distribution function of the system is made of 
a finite number N of pointwise particles. Even with a perfect 
numerical integrator, the system would necessarily undergo en¬ 
counters during which orbits feel the discreteness of the joint 
DF through its two point correlation. Note that these interactions 
need not be local but assume that the potential fluctuations are 
resonant so as to build up a secular evolution of the system. This 
effect, which is still present even in the absence of any numerical 
noise, is the effect captured by the Balescu-Lenard equation (see 
figure 1 ). 


4.1. The disc model 


The disc considered by S12 is an infinitely thin Mestel disc, for 
which the circular speed is a constant Vo independent of the 
radius. Such a model represents fairly well the observed rotation 
curve of real galaxies. The stationary background potential 
of such a disc and its associated surface density Sm are given by 


fu{R)^V^\n 


R 

A 


Ym(R) = 


^0 

IttGR ’ 


( 81 ) 


where R, is a scale parameter. Because of this scale invariance, 
the relationship between the angular momentum 7^ and the guid¬ 





Fig. 4; Surface density Z, of the tapered Mestel disc. The unit system 
has been chosen so that Vo = G = Rj = 1. Because of the tapering func¬ 
tions, the self-gravity of the disc is turned off in the inner and outer 
regions. 


ing radius Rg is straightforwardly given by 

^RgVo. (82) 

Within the epicyclic approximation, the intrinsic frequencies 
and K can be computed thanks to equations ( 22 ) and read 

y 2 

= -r ; = ^2 . (83) 

J<p 

We note that kID.^- V2, so that the Mestel disc could be seen 



Fig. 5; Contours of the initial distribution function in action-space 
(J^, Jr), within the epicyclic approximation. The contours are spaced 
linearly between 95% and 5% of the distribution function maximum. 


as an intermediate case between the Keplerian case for which 
KfQ.^-1 and the harmonic case for which /r/Q 0 = 2. The ratio 
of the intrinsic frequencies is a important parameter for the sys¬ 
tem since it will determine the location of the resonances and a 
constant ratio may introduce dynamical degeneracies. This is the 
case for the Keplerian and harmonic discs for which is a ra¬ 
tional number, as discussed below equation (69). By contrast, for 
the Mestel disc, the non-rational ratio /r/Q^ = V2 ensures that the 
potential is non-degenerate. Using the epicyclic approximation, 
the DF considered by S12 takes, as in equation (25), the form of 
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Fig. 6; Dependence of the local Q Toomre parameter with the angular 
momentun. It is scale invariant except in the inner/outer regions because 
of the presence of the tapering functions Tinner and Touter- The unit sys¬ 
tem has been chosen so that V() = G = /?, = 1. 


Fig. 7: Variations of the response matrix eigenvalues A with the WKB- 
frequency k^, for m = /wcor and two values of J^. The curve that peaks 
at large k,. is for the smaller value of 


a Schwarzschild DF, where the intrinsic frequencies are given by 
equation (83), the velocity dispersion cr^ is constant throughout 
the entire disc, and the surface density is given by St, i.e. the ac¬ 
tive surface density of the disc. Indeed, in order to accommodate 
the central singularity and the infinite extent of the Mestel disc, 
one introduces tapering functions Tinner and Touter to damp out 
the inner and outer regions, which read 

r 

Tinner(7^) = ^ 

Touter(7 0 ) — 

where v and fi control the sharpness of the two tapers and Rq 
is an additional scale parameter. The two tapers are physically 
motivated by the presence of a bulge and an outer truncation 
for the disc. Moreover, in order to reduce the susceptibility of 
the disc, we also suppose that only a fraction ^ of the disc is self- 
gravitating, with 0 < ^ < 1, so that the rest of the gravitational field 
is provided by the static halo. As a conclusion, the active surface 
density St is given by 


{RiVoY + Jl ’ 


1 - 


J(b 




(84) 


Tinner(7 cp) TQuter(7 p) , 


(85) 


where Sm is the surface density of the Mestel disc from equa¬ 
tion (81). We place ourselves in the same units system as in S12, 
so that we have Vo-G-Ri-l. The other numerical factors are 
given by cr,. = 0.284, v-A, ^-5, ^ = 0.5 and Rq- 11.5. The shape 
of the active surface density is illustrated in figure 4. The ini¬ 
tial contours of the Schwarzschild DF from equation (25) are 
shown in figure 5. For such an almost scale invariant disc, the 
local Toomre parameter, Q (Toomre 1964) 


QiJ^) = 


CTr K{ Jp) 
3.36 GS,(7^)’ 


( 86 ) 


which for Q>1 ensures the stability of the disc with respect to lo¬ 
cal axisymmetric disturbances, becomes almost independent of 
the radius, especially in the intermediate regions of the disc. As 
illustrated in figure 6, 2-1.5 between the tapers and increases 
strongly in the tapered regions. 

The expression (10) of the secular diffusion flux requires to 
sum on all the resonances /w. S12 restricted pertubations forces 
to mp-2, so that we may impose this same restriction on the 


considered azimuthal number mp. Throughout our numerical cal¬ 
culations, we will restrict ourselves to only three different res¬ 
onances which are; the inner Lindblad resonance (ILR) corre¬ 
sponding to = (-1,2), the outer Lindblad resonance 

(OLR) given by ot°^^) = (1,2) and the corotation reso¬ 
nance (COR) for which m™’^) = (0,2). Moreover, all the 

calculations in the upcoming sections have also been performed 
while taking into account the contributions from the resonances 
with ifir = +2, which were checked to be subdominant. Being able 
to perform such a restriction to the relevant resonances appear¬ 
ing in the secular flux !Ftot from equation (10) is an important 
step of the calculation. 

Returning to the fast and slow coordinates from equa¬ 
tion (14), note that the diffusion associated to the COR resonance 
amounts to diffusion along the T^-axis. Such diffusion brings 
stars from one quasi-circular orbit to another of a different ra¬ 
dius and is called radial migration. Conversely, the diffusion as¬ 
sociated to the ILR and OLR resonances exhibits a non-zero dif¬ 
fusion component in the 7r-direction. It therefore increases the 
velocity dispersion within the disc so as to heat it, while either 
decreasing (ILR) or increasing (OLR) star’s angular momentum. 


4.2. Disc amplification 

One may now study the behavior of the amplification eigenval¬ 
ues given by equation (62), thanks to which one can perform 
the improved approximation of the small denominators. For a 
given resonance m and angular momentum Jp, the amplifica¬ 
tion function kr^At^. is presented in figure 7. As equation (62) 
only depend on s^, the ILR and OLR resonances always have 
the same response matrix eigenvalues. One can also note that 
the eigenvalues A(kr) are maximum for a frequency k^aaxiJp), 
where A(kr)-A^^;i, in a region whose size is given by the width 
at half maximum Ak^. Because of the scale-invariance property 
of the Mestel disc, it is straighforward to show that Ak^ cxljjp, 
kmax “c 1 and kinf/sup i^/Jp- One can then consider the behav¬ 
ior of the amplification factor l/(l-/imax), which encodes the 
strength of the self-gravitating amplification, as shown in fig¬ 
ure 8. Note that the COR resonance is always more amplified 
than the ILR and OLR resonances, but the maximum amplifi¬ 
cation (~ 3 for the COR and ~ 1.5 for the ILR and OLR) re¬ 
mains sufficiently small, so that the susceptibility coefficients 
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Fig. 9; Map of the divergence of the total flux y'tot summed over the three resonances (ILR, COR and OLR). Red contours, for which div CFtot) < 0, 
correspond to regions from which the orhits will be depleted, whereas blue contours, for which div ('Ftot) > 0, correspond to regions where secular 
diffusion will tend to increase the value of the DF. The net fluxes involve simultaneously radial migration near (y,j,yr)~'(l-8,0), and heating near 

(y^,yo~(i,o.i). 


1/(1-Amax) 



Fig. 8; Dependence of the amplification factor 1/(1-Tmax) with the po¬ 
sition in the disc. The amplification associated to the COR is always 
larger than the one associated to the ILR and OLR. 


from equation (4) are not dominated only by the self-gravitating 
amplification. 

4.3. Computing the diffusion fiux 

Given the knowledge of the eigenvalues, it is now straight¬ 
forward to compute the susceptibility coefficients within the 
improved approximation of the small denominators thanks to 
equation (77), where the integration on kr is performed for 
[kinf; ksup]. One can then compute the associated drift 
and diffusion coefficients respectively given by equations (79) 
and (80). The diffusion flux y^tot defined in equation (10) im¬ 
mediately follows, where the sum on m is restricted only to the 
three resonances ILR, COR and OLR. In Appendices E and F, 
we discuss two specific properties of such a truncated Mestel 
disc, namely the cancellation between the radial components of 
the diffusion and drift elements (the Schwarzshild conspiracy, 
Appendix E) and the natural and intrinsic presence of a tempo- 



Fig. 10; Overlay of the WKB predictions for the divergence of the 
diffusion flux div(!Ftot) and the differences between the initial and the 
evolved DF in S12’s simulation. The opaque contours correspond to 
the differences in the action-space for the DF in S12 between the time 
tsi2=1000 and rsi 2 = 0 (see the upper panel of S12’s figure 10). The 
red opaque contours correspond to negative differences, so that these re¬ 
gions are emptied from their orbits, whereas blue opaque contours cor¬ 
respond to positive differences, i.e regions where the DF has increased 
through diffusion. The transparent contours correspond to the predicted 
values of div ('Ftot) using the same conventions as in figure 9. Note the 
overlap of the predicted transparent red and blue contours with the mea¬ 
sured solid ones. 


ral frequency bias (Appendix F), which both enlighten the subtle 
arbitrations between the different resonances. 

Finally let us compute the divergence of this flux, div (’Ftot) 
given by equation (11) in order to compare quantitatively the 
WKB predictions with the results from S12’s simulations. Fig¬ 
ure 9 represents the contours of div (Ftot)- A comparison of our 
WKB predictions with the results from S12’s simulation are il¬ 
lustrated in the figures 10 and 11. In figure 9, red contours 
correspond to regions for which div (y'tot)<0, so that thanks 
to equation (11) they are associated to action-space regions 
where the WKB Balescu-Lenard equation predicts a decrement 
of the DF during secular evolution. In contrast, blue contours 
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Fig. 11: Overlay of the WKB predictions for the divergence of the dif¬ 
fusion flux div (!Ftot) on top of the contours of the DF in action-space 
measured in the S12 simulation. The black background contours are the 
levels contours of the DF at time fsi 2 = 1400 (see the lower panel of fig¬ 
ure 7 of S12). These contours are spaced linearly between 95% and 5% 
of the DF maximum and exhibit clearly the appearance of a resonant 
ridge. The colored transparent contours correspond to the predicted val¬ 
ues of div CFtot) using the same conventions as in figure 9. Note that the 
developed late time ridge is consistent with the predicted depletion (red) 
and enrichment (blue) of orbits. 


are associated to regions for which div (ytot) >0, so that the 
DF will increase there. The overall picture involves two com¬ 
peting processes: i) the beginning of a ridge forming towards 
and ii) the formation of an over density near 
(J^, 0). Point i) is in fact consistent with both the early 

time measurement of S12 as shown in hgure 10 and the late time 
measurement of S12 as shown in hgure 11. These qualitative 
agreements are in fact surprisingly good, given that the WKB 
theory is approximate and was only estimated for f = 0^. Interest¬ 
ingly, the early time measurement from hgure 10 also displays 
a hint of an over-density on the Jr-Q axis, in agreement with 
point ii), while the late time measurement suggests that the over 
density has split, with a hint of a second ridge forming. The time 
evolution of equation (2) is likely to explain why this over den¬ 
sity on the axis seems to split, and why the ridge gets amplihed. 

From hgure 9, we explicitly compute div (’Ftot), so that we 
may now study the typical timescale of collisional relaxation pre¬ 
dicted by the WKB Balescu-Lenard equation. This is the purpose 
of the next section. 

4.4. Physical timescales 

Given estimates of the diffusion Hux, one can explicitly compute 
the collisional timescale, i.e. the timescale for which the hnite- 
N effects become signihcant. Indeed, the larger the number N of 
particles, the later these effects will come into play. One should 
note that our writing of the Balescu-Lenard equation (2) is inde¬ 
pendent of the number N of particles. However, when correctly 
dimensionalized, this kinetic equation takes the form 

^+L[F] = icBL[F], (87) 

where L is the operator of pure advection, and Cbl is the Balescu- 
Lenard collisional operator, i.e. the r.h.s of equation (2). Equa¬ 
tion (87) underlines the fact that the collisional term is associ¬ 
ated to a kinetic Taylor expansion in the parameter e-\ jN 


Within the angle-actions coordinates, the advection operator is 
immediately given by 

L = . (88) 

oO 

Because we have assumed that F is always quasi-stationary, so 
that F - F(J, t) (adiabatic approximation), one has L[F'] = 0. 
We now introduce the time 



(89) 


so that equation (87) immediately becomes 

^ = CblIF] . (90) 

Equation (90) corresponds to a rewriting of the Balescu-Lenard 
equation, where N is not present anymore. This will allow us to 
quantitatively compare the time during which the S12 simulation 
was run to the diffusion time predicted by our WKB Balescu- 
Lenard formalism. In order to ease this comparison, we place 
ourselves in the same units system as the one used by S12. Eig- 
ure 7 of S12 for which the ridge was observed was obtained with 
the parameters N-50x 10® and Afsi 2 = 1400. Using the rescaled 
time introduced in equation (89), one obtains that S12 observed 
the resonant ridge after a time Atsi 2 = Afsi2/A^-3xl0 One 
can then compare this time, with the typical time required to ob¬ 
tain a resonant ridge within our WKB formalism. Given the map 
of div(y'tot) described in section 4.3, one can estimate the typi¬ 
cal time for which this flux could lead to the features observed 
in SI2. The contours presented in the hgure 7 of S12 are sepa¬ 
rated by a value of 0.1 xF™”, where ^ 0.12 corresponds 
to the maximum of the normalized DE from equation (25). As a 
consequence, to observe the resonant ridge, the DE should typ¬ 
ically change by a value of the order of AFo-0.1xF™“. Erom 
hgure 9, one can note that the maximum value of the divergence 
of the Hux is given by |div (’Ftot) I - 0.06. Einally, thanks to equa¬ 
tion (11), one may write the relation AFq-Atwkb |div(ytot)l, 
where Atwkb is the minimal time during which the WKB 
Balescu-Lenard equation has to be considered in order to de¬ 
velop a ridge. Thanks to the previous typical numerical values, 
one obtains that ArwKB-3xl0 '. When comparing these two 
typical times, Atsi 2 the duration during which S12 simulation 
was performed and Atwkb the duration required to observe secu¬ 
lar diffusion in the WKB Balescu-Lenard formalism, one obtains 
the order of magnitude 


Atsi2 

Atwkb 


10 ^ 


(91) 


Hence the direct application of the WKB-limited Balescu- 
Lenard equation does not allow us to predict the observed 
timescale for the diffusion features in simulations. Indeed, the 
timescale of collisional diffusion predicted by this WKB formal¬ 
ism seems much larger than the time during which the numerical 
simulation was performed. This discrepancy is also strengthened 
by the use of a softening length in numerical simulations, which 
induces an effective thickening of the disc, so as to slow down the 
collisional relaxation. A possible explanation for this timescale 
discrepancy is discussed in the next section. 


4.5. Interpretation 

In order to interpret S12 simulation under the light of a colli¬ 
sional secular diffusion equation, such as the Balescu-Lenard 
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equation and its WKB limit, one should first note the undisputed 
presence of collisional effects in S12’s simulation. Indeed, figure 
2 of S12 shows that when the number of particles of the simula¬ 
tion is increased, the strength of the density fluctuations are de¬ 
layed, which in turn is likely to be related to the amplitude of the 
secular features. The larger the number of particles, the later the 
effect of secular diffusion. Such dependence illustrates the fact 
that discreteness effects do play a role in the secular diffusion 
observed in S12. 

Sellwood & Kahn (1991) have argued that a sequence of 
causally connected swing amplified transients could occur sub¬ 
ject to a (possibly non local) resonant condition between succes¬ 
sive spirals waves. The Balescu-Lenard formalism captures pre¬ 
cisely such sequences - in as much as it integrates over dressed 
correlated potential fluctuations subject to relative resonant con¬ 
ditions, but does not preserve causality nor resolve them on dy¬ 
namical timescales. The exact initial phases are not relevant in 
the Balescu-Lenard formalism: see Appendix A for a sketch of a 
full derivation which makes this point clear. 

The timescale discrepancy observed in equation (9 1 ) might 
be driven by the incompleteness of the WKB basis. Indeed, equa¬ 
tion (26)’s basis - thanks to which the susceptibility coefficients 
were evaluated in equation (77) - does not form a complete set, 
as it can only represent correctly tightly wound spirals. It also en¬ 
forces local resonances, and does not allow for remote orbits to 
resonate, or wave packets to propagate between such non local 
resonances. The seminal works from Goldreich & Lynden-Bell 
(1965); Julian & Toomre (1966); Toomre (1981) showed that 
any leading spiral wave during its unwinding to a trailing wave 
undergoes a significant amplification, coined swing amplifica¬ 
tion. Because it involves open spirals this mechanism is not cap¬ 
tured by the WKB formalism. This additional amplification is 
expected to increase the susceptibility of the disc and therefore 
accelerate secular diffusion (both drift and diffusion), so that the 
timescales discrepancy from equation (91) will become less re¬ 
strictive. Following the notations from Toomre (1981), the trun¬ 
cated Mestel disc considered in the S12 simulation corresponds 
to 2^ 1.5 and X = 2, so that figure 7 from Toomre (1981) shows 
that significant swing amplification (of order ~ 10) may be ex¬ 
pected. It has also been claimed (Toomre & Kalnajs 1991) that 
swing amplified shot noise in the shearing sheet approximation 
would behave like significantly heavier macro-particles. Such an 
amplification would keep a dependence of the secular response 
with the total number N of particles, but would reduce signifi¬ 
cantly the effective number of particles. 

The Balescu-Lenard WKB limit seems to capture qualita¬ 
tively the main features of the initial diffusion process in action 
space (as discussed in Section 4.3), but falls short in predicting 
the relevant timescale. The remaining questions are therefore: 
what is the exact impact of swing amplification? Can it explain 
the timescale discrepancy? 

5. Conclusion 

We implemented the inhomogeneous Balescu-Lenard equa¬ 
tion (2) for an infinitely thin galactic disc using two main approx¬ 
imations. We first assumed the disc to be tepid. We could then 
use the epicyclic approximation which allowed for an explicit 
mapping between the physical coordinates (x, v) and the angle- 
actions coordinates (6, J) via equation (24). Our second approx¬ 
imation relied on the introduction of the tightly wound basis el¬ 
ements from equation (26). Because of the corresponding WKB 
approximation, we obtained in equation (62) a diagonal response 
matrix, so that gravity is effectively treated locally. The associ¬ 


ated scale-decoupling hypothesis yields a crucial restriction to 
only local resonances, as shown in equation (72). We then de¬ 
rived in equation (77) a simple quadrature for the susceptibility 
coefficients, given by equations (B.7) and (B.8) for the bare ones. 
Thanks to this restriction to local resonances, we were also able 
to write the drift and diffusion coefficients as simple quadratures 
in equations (79) and (80). 

These simple expressions derived within the WKB formal¬ 
ism yield, to our knowledge, a first non trivial explicit expression 
for the Balescu-Lenard diffusion and drift coefficients. They are 
certainly useful to provide insight into the physical processes 
at work during the secular diffusion of a self-gravitating discrete 
disc. Moreover, modulo the restriction to the three physically mo¬ 
tivated resonances ILR, COR and OLR, our WKB formalism can 
be used for quantitative comparisons to numerical experiments 
such as the one presented in Sellwood (2012). It considered a sta¬ 
ble isolated Mestel disc sampled by pointwise particles, whose 
secular evolution is induced by finite-A effect ideally captured 
by the Balescu-Lenard equation. 

The straightforward calculation in the WKB limit of the di¬ 
vergence of the full diffusion flux, div ('Ftot), (illustrated in fig¬ 
ures 9, 10 and 11), recovered most of the secular features ob¬ 
served in SI2. This qualitative agreement is impressive, given 
the level of approximation involved in the WKB limit. The hints 
for the formation of a ridge - depletion and enrichment of or¬ 
bits along a preferred direction - is qualitatively consistent with 
the findings of S12 and Fouvry & Pichon (2015); Fouvry et al. 
(2015), without postulating additional assumptions about the 
source of fluctuations^. 

The comparison of the collisional time predicted in the WKB 
limit (equation (91)) to the diffusion time of S12 simulation, 
highlights nonetheless a significant quantitative overestimation. 
We provided a possible explanation which relies on the intrin¬ 
sic limitations of the WKB formalism, as it cannot account for 
swing amplification, during which unwinding transient spirals 
are strongly amplified. This additional amplification, which in¬ 
volves explicitly non local wave absorption and emission, may 
be the missing contribution required to reconcile quantitavely 
our predictions and the simulation. One venue will be to com¬ 
pute numerically exactly equations (8) and (9) in action space - 
without assuming tightly wound spirals or epicyclic orbits - with 
a complete basis This is the topic of an upcoming numerical 
investigation (Fouvry et al. 2015). 

Should this complementary investigation explain the 
timescale mismatch, one would be in a stronger position to 
validate the accuracy of A-body schemes to correctly capture 
secular evolution of discrete self-gravitating cold discs over very 
long timescales. This would clearly be a worthy assessment of 
such schemes relying on the Balescu-Lenard theory. Once the 
above described conundrum is resolved, we also will be able 
to evolve over secular times the Balescu-Lenard equation and 
predict the full cosmic time evolution of such discrete discs. 
This may also contribute to solving the timescale discrepancy. 

In closing, beyond the application presented in section 4, the 
above developed tightly wound Balescu-Lenard formalism may 
for instance describe the secular diffusion of giant molecular 
clouds in galactic discs (which in turn could play a role in mi¬ 
gration driven metallicity gradients and disc thickening), the sec- 

^ In contrast, the formalism of secular forcing presented 
in Fouvry & Pichon (2015); Fouvry et al. (2015) postulated a partially 
ad-hoc shape of the perturbation power spectrum, see equation (F.2). 

^ An alternative middle ground would be to account for non local in¬ 
terferences of the WKB wave packets, given by equation (26), which in 
turn would allow for non local resonances to come into play. 
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ular migration of planetesimals in partially self-gravitating proto¬ 
planetary discs, or even the long-term evolution of population of 
stars, gas blobs and debris near the Galactic center. Such topics 
will be subject to further investigations. 

Acknowledgements. JBF thank the Institute of Astronomy, Cambridge, for hos¬ 
pitality while this investigation was completed. JBF and CP also thank the the¬ 
oretical physics sub-department, Oxford, for hospitality and the CNRS-Oxford 
exchange program for funding. We thank Donald Lynden-Bell, James Binney, Si¬ 
mon Prunet, Walter Dehnen, John Magorrian and Mir Abbas Jalali for their feed¬ 
back. This work is partially supported by the Spin(e) grants ANR-13-BS05-0005 
of the French Agence Nationale de la Recherche (http: //cosmicorigin. org) 
and by the LAB EX Institut Lagrange de Paris (under reference ANR-IO-LABX- 
63) which is funded by ANR-1 l-IDEX-0004-02. 


References 

Balescu, R. 1960, Physics of Fluids, 3, 52 

Binney, J. & Tremaine, S. 2008, Galactic Dynamics: (Second Edition), Princeton 
Series in Astrophysics (Princeton University Press) 

Bom, M. 1960, The Mechanics of the Atom (F. Ungar Pub. Co.) 

Chandrasekhar, S. 1942, Principles of Stellar Dynamics (University of Chicago 
Press) 

Chandrasekhar, S. 1949, Rev. Mod. Phys., 21, 383 

Chavanis, P.-H. 2007, Physica A Statistical Mechanics and its Applications, 377, 
469 

Chavanis, P.-H. 2012a, Physica A Statistical Mechanics and its Applications, 391, 
3680 

Chavanis, P.-H. 2012b, Physica A Statistical Mechanics and its Applications, 
391, 3657 

Chavanis, P.-H. 2012c, European Physical Journal Plus, 127, 19 
Chavanis, P.-H. 2013a, European Physical Journal Plus, 128, 126 
Chavanis, P.-H. 2013b, A&A, 556, A93 

Chavanis, P.-H. & Lemou, M. 2007, European Physical Journal B, 59, 217 

Daubechies, I. 1990, Information Theory, IEEE Transactions on, 36, 961 

Earn, D. J. D. & Lynden-Bell, D. 1996, MNRAS, 278, 395 

Fouvry, J.-B., Binney, J., & Pichon, C. 2015, submitted 

Fouvry, J.-B. & Pichon, C. 2015, MNRAS, 449, 1982 

Fouvry, J.-B., Pichon, C., & Prunet, S. 2015, MNRAS, 449, 1967 

Fouvry, J. B. et al. 2015, in prep 

Gabor, D. 1946, Electrical Engineers - Part III: Radio and Communication Engi¬ 
neering, Journal of the Institution of, 93, 429 
Goldreich, P. & Lynden-Bell, D. 1965, MNRAS, 130, 125 
Goldstein, H. 1950, Classical mechanics (Addison-Wesley) 

Gradshteyn, 1. S. & Ryzhik, 1. M. 2007, Table of integrals, series, and products 
(Elsevier/Academic Press, Amsterdam) 

Heyvaerts, J. 2010, MNRAS, 407, 355 

Jeans, J. 1915, Monthly Notices of the Royal Astronomical Society, 76, 70 
Jeans, J. 1929, Astronomy and Cosmogony (Cambridge Univ. Press) 

Julian, W. H. & Toomre, A. 1966, ApJ, 146, 810 
Kalnajs, A. J. 1965, Ph.D. thesis (Harvard University) 

Kalnajs, A. J. 1976, ApJ, 205, 745 

Klimontovich, I. 1967, The statistical theory of non-equilibrium processes in a 
plasma (M.I.T. Press) 

Landau, L. 1936, Phys. Z. Sowj. Union, 10, 154 
Lenai'd, A. 1960, Annals of Physics, 10, 390 

Lin, C. C. & Shu, F. H. 1966, Proceedings of the National Academy of Science, 
55,229 

Liouville, J. 1837, J. Math. Pures AppL, 1, 16 
Lynden-Bell, D. 1967, MNRAS, 136, 101 
Lynden-Bell, D. 1979, MNRAS, 187, 101 
Lynden-Bell, D. & Kalnajs, A. J. 1972, MNRAS, 157, 1 
Mestel, L. 1963, MNRAS, 126, 553 

Palmer, P. 1994, Stability of Collisionless Stellar Systems: Mechanisms for the 
Dynamical Structure of Galaxies, Astrophysics and Space Science Library 
(Springer Netherlands) 

Palmer, P. L., Papaloizou, J., & Allen, A. J. 1989, MNRAS, 238, 1281 
Pichon, C. 1994, Ph.D. thesis (University of Cambridge) 

Pichon, C. & Cannon, R. C. 1997, MNRAS, 291, 616 

Rosenbluth, M., MacDonald, W., & Judd, D. 1957, Phys. Rev., 107, 1 

Sellwood, J. A. 2012, ApJ, 751, 44 

Sellwood, J. A. & Kahn, F. D. 1991, MNRAS, 250, 278 

Toomre, A. 1964, ApJ, 139, 1217 

Toomre, A. 1977, ARA&A, 15, 437 

Toomre, A. 1981, in Structure and Evolution of Normal Galaxies, ed. S. M. Fall 
& D. Lynden-Bell, 111-136 

Toomre, A. & Kalnajs, A. J. 1991, in Dynamics of Disc Galaxies, ed. B. Sun- 
delius, 341 

Vlasov, A. 1938, Zh. Eksp. i Teor Fiz., 8, 291 
Weinberg, M. D. 1993, The Astrophysical Journal, 410, 543 


Appendix A: Sketch of Balescu-Lenard Derivation 

Two derivations of the inhomogeneous Balescu-Lenard equation 
have been presented in the literature. The first one (Heyvaerts 
2010) is based on the appropriate truncation at the order 1/N 
of the BBGKY hierarchy. The second (Chavanis 2012a) relies 
on the Klimontovich equation, using a quasilinear approxima¬ 
tion. We now briefly sketch the derivation presented in Chavanis 
(2012a). We consider an isolated system of N particles in inter¬ 
action, of mass OT= 1, in a physical space of dimension d. Their 
dynamics is entirely described by Hamilton’s equations 

djc, dH di), dH tA li 

df dvi ’ df dxi ’ 

where the Hamiltonian of the system is given by 

N J 

^ -H ^ u(\Xi-Xj\) . (A.2) 

1=1 i<j 

Here u(\Xi-Xj\) is the binary potential of interaction. In the grav¬ 
itational case, it satisfies u{x)--Gl\x\. One can now introduce 
the discrete distribution function Fi{x, v, t) defined as 


Fi(x, D, f) = ^ 5 d(a:-J(:,(0) 5d(«-«;(0) , 
1=1 

along with the corresponding potential 
ij/i^x, t) - J'dx'dv' u{\x-x'\) T'd(jc', t). 


(A.3) 


(A.4) 


One can show that F^ satisfies the Klimontovich equa¬ 
tion (Klimontovich 1967) 


dFj ^ dHj dFj dHd dFj 
dt dv dx dx dv 
where we have defined the Hamiltonian as 

1 9 

77d = 2 " t). 


(A.5) 


(A.6) 


At this stage, it is important to note that the Klimontovich 
equation (A.5) contains exactly the same information as the 
Hamilton equation (A.l). We now introduce the smooth dis¬ 
tribution function F(x,v,t)-{Fd{x,v,t)), corresponding to an 
average of F^ over a large number of initial conditions. One 
can then write F^-F+dF, where 6F denotes fluctuations 
around the smooth distribution. In a similar way, we introduce 
ilf(x, V, t) - {xjJi{x, V, t)), so that - i/z+di/z. We have therefore de¬ 
composed the discrete distribution function F^ into a smooth 
component F that evolves slowly with time, whereas the fluc¬ 
tuating component 6F evolves more rapidly. As a consequence, 
when considering the evolution of the fluctuations, one can as¬ 
sume the smooth distribution to ht frozen. Using this timescale¬ 
decoupling approach, one can use the angle-actions coordinates 
{0\,J\) associated with the quasi-stationary smooth potential ij/ 
to describe the fast evolution of the fluctuations. Using these 
decompositions and this change of coordinates, equation (A.5) 
takes the form of two evolution equations 


d6F ^ d6F ddi// dF 

----- 

dt dOi 56ii 57i 

and 


— 

dt 


d 

dJx 


6F 


ddijj 


(A.7) 


(A.8) 
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where we have introduced the intrinsic frequencies of the system 
Qi as in equation (1), and where {.) denotes an angle 

average. Because of our timescale-decoupling approach, we may 
neglect the time variation of F(Ji, t) in the calculation of the col¬ 
lision term (adiabatic approximation). It implies that F may be 
treated as a constant in equation (A. 7), because F evolves on a 
(relaxation) timescale much larger than the (dynamical) time cor¬ 
responding to the evolution of 6F. In order to be valid, this ap¬ 
proximation requires to have A » 1. Finally, we also assume that 
the distribution F remains Vlasov stable, so that its evolution is 
only governed by correlations and not by dynamical instabilities. 

The first step of the derivation of the Balescu-Lenard equa¬ 
tion is then to study the short timescale evolution equation (A. 7). 
One defines the Fourier-Laplace transform of the fluctuation 6F 
as 


5F„„(/i,wi) = J ^dfc-'''"> ®>-“''>dF(0i,/i,f), (A.9) 

valid for Im(a;i) sufficiently large. Similarly to equation (6), one 
also defines the spatial Fourier tranform of the initial value as 

^e-''">®>(5F(0i,7i,O). (A.IO) 

Thanks to these transformations, equation (A. 7) may be rewrit¬ 
ten under the form 


SFmiiJ 


niydFIdJy ~ 5F'mj(7i,0) 

---Wi)+T^-Pi-r 


(A. 11) 


We now use the basis elements introduced in equation (3), so that 
we may decompose the potential fluctuations under the form 

6ilr(0uJi,t) = 2flp(f)t^''’^(0i,7i). (A. 12) 

p 


We introduce the Laplace transform of ap(t) as 


Using equation (A. 1 1), one immediately obtains that 

mi-dFIdJi .r 

-Pi-(/1 , Wl) U UOJl)} 

{6F,„^{J 1,0) 6iJ/m2iJ 1 ,(^ 2 )) 

/(mi'Qi-wi) 

In equation (A. 17), the first term corresponds to the self¬ 
correlation of the potential whereas the second term corresponds 
to the correlations between the potential fluctuations and the dis¬ 
tribution function at time t-Q. Each of these terms must then be 
considered one at a time. Assuming that there is no correlation in 
the initial phases, one can show (see Appendix C from Chavanis 
2012a) that 

<dA«.(/i, 0) 6F,„,(J2, 0)) = SMi-Ji) F(Ji). (A. 18) 


Hence, given equation (A. 1 5), one can rewrite the first term of 
equation (A. 17) under the form 




1 




d/3 


FiJi) 


1 






,iJl,J3,OJ2) (/M3-Q3-Wi)(/M3-Q3-HW2) 


3(/i,/ 3, Wl) 


. (A.19) 


If we consider only the contributions that do not decay in time, 
one can perform the substitution 


1 

(ni3-Q3-Wi)(ni3-Q3H-W2) 


^i27T)^6DioJi+a>2) Soimi-Sli-oJi). 


Starting from equation (A.19), thanks to the previous substitu¬ 
tion and using the fact that .m, (/, / 3 , -wi) = .m, (/ 1 , / 3 , wi )*, 

one can show that the first contribution from equation (A. 16) 
takes the form 


-1 


ap(o)i) - dtap(t) e 


iciJit 


(A. 13) 


(m2-il2-0Ji) 




27r 

m\,m2 




d /2 mi 


itiydFIdJi 


Let us then take the inverse Fourier transform of equation (A. 1 1), 
multiply by \) and integrate over 0i and J\ (using the 

property that djc dt; = d^i d/i). One gets 

5^(cui) = -(27r)‘5)/-M(wi)]^;2jd/2T 


ffly • ^2l — (jJ\ 
5D{m2-Sl2-OJi) 


i,m2 (/l,/2, Wl)h 


FU 2 ) . (A.20) 


<JFm,(/2,0) , 

(J2) , 


The last step of the calculation is to use the Landau prescription 
a)^a>+iQ^ along with Plemelj formula 


1 


(A. 14) Jc + i0+ 




induix), 


(A.21) 


where the response matrix M is given by equation (5). Equa¬ 
tion (A. 14) can be rewritten under the form 


(/i,wi) = -(27r)' 




1 


5F„,(/2,0) 


iDm2,m2iJl,J2,OJ\) i{m2-^2-<X>\) 
(A.15) 


where the susceptibility coefficients have been introduced in 
equation (4). One can now compute the collision term appear¬ 
ing in the r.h.s of equation (A.8). It requires us to evaluate 

lsF^] = y [^^im2 

\ ^^ 1 / 2^ 27r 

X <dFm,(/l,Wl)#mj(/l,W2)). (A.16) 


where V is the principal value. One can then rewrite equa¬ 
tion (A.20) under the form 


SF^TTiln) 
dOi /i 


m\,m2 


d/2 mi 


5u{mi-^\-m2-il2) 

,m2 (/i,/2,mi-fli) 


dF\ 

xK— 


(A.22) 


Similarly, one can rewrite the second term of equation (A. 1 7) 
under the form 


<^Fmi (/1,0) (/1. tL>2) ) 

/(mi'Qi-wi) 


,2^D(^i+^2)<5D(mrni-wi) 

“( 271 ) -—-—— --- r(Ji). 


'D„ 




(A.23) 
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Using symmetries of the matrix starting from equa¬ 

tion (A.23), one can show that the second contribution from 
equation (A. 16) finally takes the form 



-7T (In) 


r iJim 

m\,m2 


,m2 


X 



FUx)- 


(A.24) 


Combining the two contributions obtained in equations (A.22) 
and (A.24), one can rewrite equation (A. 16) under the form 


6F 


86x1/ 


-n {2n) 


mi ,m2 




6u(ml■^ll-m2■^l2) 
^ju(/i,f)U(/2,0. (A.25) 


Hence, using the slow evolution equation (A. 8), we recover the 
Balescu-Lenard equation introduced in equation (2). As a final 
remark, one must note that on the short dynamical timescale, the 
evolution is governed by equation (A.7), which involves the fluc¬ 
tuating components 6F and 6x1/ of the distribution function and 
the potential. In contrast, on the long secular timescale, the evolu¬ 
tion is governed by equation (A. 8), which after an angle-average 
only involves the mean distribution function F. Indeed, thanks 
to the ensemble average, all the cross-correlations between the 
fluctuations 6F and 6x1/, as in equations (A. 1 8), (A. 1 9) and (A.23) 
can be expressed in terms of the underlying smooth distribution 
function F only, so that the fluctuating components are absent 
from the secular Balescu-Lenard collision operator from equa¬ 
tion (A.25). 


susceptibility coefficients from equation (4), one must necessar¬ 
ily rely on Kalnajs’ matrix method (Kalnajs 1976). It is however 
possible to express using the potential basis. Indeed, for a 

fixed value of X 2 , we consider the function Xi u(xi-X 2 )- One 
can then decompose this function on the basis elements ^^p\xi), 
so that we may write 

u(xi-X 2 ) = ^ ap(x 2 ) xl/^'’\xi) , (B.3) 

p 


where it is important to note that the basis coefficients ap{x 2 ) are 
functions of JC 2 . Thanks to the biorthogonality property of the 
basis detailed in equation (3), one can obtain the expression of 
the coefficients ap{x 2 ) which reads 


ap{x2) = -J' 1 a:i u(xi-X2)p^^^*(xi) 

-- Jdxi u(xi-X 2 )p^^\xi) 

^-xp(P>{r2), 


(B.4) 


where we used the fact that m is a real function. Hence we Anally 
obtain 

u{Xi-X2) = - ^Xp^PHxx)xp^’’^*{X2) ■ (B.5) 

p 


Taking appropriately a Fourier transform with respect to the an¬ 
gles as in equation (B.2), we obtain 

A„,.„,(7i,/2) (B.6) 

p 


Appendix B: Turning off coiiective effects 

The reference Chavanis (2013b) (see also Appendix A of 
Chavanis (2012a)) considers the inhomogeneous Balescu- 
Lenard equation without collective effects. This collisional 
kinetic equation is the equivalent of the Landau equation for 
inhomogeneous systems. It can be straightforwardly obtained as 
an approximation of the Balescu-Lenard equation (2). Indeed, 
one has to make the substitution from the dressed susceptibil¬ 
ity coefficients 11Dm,,m 2 (Ji, J 2 , <F) to the bare ones given by 
Am,.m 2 iJi 2 J 2 ), so that the inhomogeneous Balescu-Lenard equa¬ 
tion without collective effects {i.e. the inhomogeneous Landau 
equation) is given by 


of 




d/2 doimi ■ Qi -m2■ Q2) 


d d 

X |Affl,,„,(/i,/2)|^(mi- —-m2-— |F(/i,f)F(/2,f) 


, (B.l) 


where the coefficients are associated to the Fourier 

transform in angles of the interaction potential (Pichon 1994; 
Chavanis 2013b) and read 


^m\ ,m2 (/l,/2) = 

^J i0id02 u{\x{0x,Ji)-x{02, / 2 )|)e-'('«> «‘-”7-«4 ^ (b. 2 ) 

where u(x) is the binary interaction given by u{x) = —Gl\x\ in the 
gravitational case. A key remark at this stage is that the expres¬ 
sion (B.2) does not require to introduce biorthogonal basis ele¬ 
ments as in equation (3), whereas in order to estimate the dressed 


This fairly simple relation allows us to express the bare suscepti¬ 
bility coefficients A^,using the potential basis. One does not 
need anymore to perform a Fourier transform in angles of the 
interaction potential, because the resolution of Poisson’s equa¬ 
tion has been implicitly hidden in the effective construction of 
the basis elements. Neglecting the collective effects in the ex¬ 
pression (4) of the dressed susceptibility coefficients amounts to 
taking M(w) = 0, so that we obtain 


1 

,/«2 


~ ^m\,m2^JJ2^ ■ 

w/ocoll. 


(B.7) 


The negative sign in equation (B.7) plays no significant role in 
the kinetic equations, since one has to make the substitution of 
the square modulus l/|/)|^i-^ |Ap in the Balescu-Lenard equa¬ 
tion (2), to obtain the inhomogeneous Landau equation (B.l). Us¬ 
ing our WKB basis, one can proceed similarly as in the dressed 
case, by limiting oneself only to local resonances. Equation (76) 
therefore becomes in the bare case 


Ami,mi jj) 


Itt RI 


X 


f 

Ji/i 


dkf ^ ^ 

l/o-i 




(B.8) 


This expression of the bare susceptibility coefficients can be 
straightforwardly obtained from the dressed ones by imposing 

^ This also explains why the factorization assumption 
A„,,,m,(/i,/ 2 ) = A,„,(/,)A„, 2 (/ 2 ) used in Chavanis (2007) does 
not hold. 
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/1 = 0. One should note that because of the absence of the amplifi¬ 
cation term the approximation of the small denomina¬ 

tors cannot be used. Hence a physically motivated evaluation of 
the expression (B.R) becomes more subtle to perform, especially 
because of the possibly important role that the Coulomb loga¬ 
rithm 1 jkr might play. Finally, one can note that even for exactly 
local resonances, Le nii - m 2 and R\ -R 2 , the use of our WKB 
formalism allowed us to obtain non-diverging bare susceptibility 
coefficients. On the other hand, one could try to estimate the bare 
susceptibility coefficients starting from equation (B.2), i.e. with¬ 
out using any potential basis. Using the polar coordinates {R, (f>), 


one has to compute 



a„,,™,(7i, 72) = - ^ Jde^defde^de^ 



V - 


(B.9) 

^R\ +R\-2RiR2 cos(0i -(^ 2 ) 

By azimuthal symmetry, it is 
show (Pichon & Cannon 1997) that 

straightforward to 

^m\ ,m2 0 . 


(B.IO) 


Hence the different m^-modes of the Ami,m 2 coefficients are in¬ 
dependent. This result is identical to what was obtained in equa¬ 
tion (65) for the dressed case. In order to illustrate the regular¬ 
izing role of the WKB basis from equation (26), we will place 
ourselves in the context of an extremely tepid disc, and therefore 
assume j} -J^-Q. Thanks to the epicyclic mapping from equa¬ 
tion (24), one can drop all the dependences on Or appearing in R 
and (f). Equation (B.9) then immediately implies that 


of the bare susceptibility coefficients derived from the Fourier 
transform of the interaction potential as in equation (B.2), when 
restricted to exactly local resonances becomes logarithmically di¬ 
vergent. This divergence, which is observed in the case of local 
resonances, is induced by the interaction of singular orbits. It is 
important to note that this divergence was not observed in equa¬ 
tion (B.8) when computing the bare susceptibility coefficients 
using the WKB basis from equation (26). This implies that the 
WKB basis is not complete. For a complete biorthogonal basis, 
the expression (B.6) of the A,„^ m 2 coefficients is an exact expres¬ 
sion. However, our calculation shows that the scale-decoupled 
WKB basis we used, possesses a subtle regularizing incomplete¬ 
ness which allows to get rid of the diverging contributions to the 
coefficients in the limit of exactly local resonances. 

Appendix C: The Schwarzschild DF case 

When considering a Schwarzchild distribution function as in 
equation (25), while relying on the approximation of the small 
denominators from equation (78), one can explicitly perform 
the remaining integration on the radial action jj in the expres¬ 
sions (79) and (80) of the drift and diffusion coefficients. We now 
detail this explicit calculation. For such a Schwarzschild distribu¬ 
tion function, it is straightforward to check that the gradients of 
the distribution function with respect to the actions are given by 

dF K dF d 

_ _ JT . _ JT 

dJr cr^ ’ dJip dJ^ 

Using the expression of the susceptibility coefficients from equa¬ 
tion (78), after some simple algebra, one can rewrite the drift 
coefficients from equation (79) under the form 




(C.l) 


m[ — m 2 — 0 , 

and the bare susceptibility coefficients are given by 

G 


,iRuQ,R2,0) ^ 


X<;f< 




^Jr]-\-RI- 2 RiR 2 cos(fl^- 


(B.ll) 


(B.12) 


In this illustrative limit, we have by construction no contributions 
from the ILR and OLR resonances for which mr + 0. After an 
immediate change of variables, using becomes 


A/hj {.J 1) — 


d 




Inl 


QE 


dJ, exp 
d 


K - 

- — ■’r 
cri 


J 


— 


-Jt 


dJ, 


1 2 
cri 


(C.2) 


Similarly, the diffusion coefficients from equation (80) take the 
form 


Dmi(Jl) — 



a-t 


J 



(C.3) 


In equations (C.2) and (C.3), in order to shorten the notations, 
we introduced the function J}) defined as 


^Q,m^,0,mi/,(Fu 0^ R2 t 0} — 


The resonance condition from equation (68) takes the form 
m^Q.^{Ri)-m^Q.^{R 2 ), because we restricted ourselves in equa¬ 
tion (B. 1 1) to the case = 0. Hence assuming that the function 
Rg i-> Q.^(Rg) is a monotonic function, one has to satisfy the con¬ 
straint Ri -R 2 , so that we are restricting ourselves only to lo¬ 
cal resonances as in equation (72). For such resonances, equa¬ 
tion (B. 1 3) becomes 


G r 

— dA- 


f'm^A 


sJ(Ri-R2)^+2RiR2(l-cos(A)) 


(B.13) 


Ri,0) - - 


G r 

-^ dA 


Vl-cos(A) 


(B.14) 


gmi(J,fi, Jr) — 


1 QE (Ak^)^ 

1 


1 ” ■^max 


l2 


Jr 


where we used the same shortened notation for the resonant fac¬ 
tor \lim\-il\)' as in equation (E.3). In addition to the integration 
formula (57), we may also rely on the additional identity 


X-FOO 

dJrJre-<^^'jl[lis[J^ = 


e 




P 

- — + \ + \mr 

2a 


J- rr^ 


£ 

2a 


ti 


2a 


(C.4) 


where a > 0, y8 > 0, and 6 Z. In analogy with the dehnition from 

equation (58), we also introduce max as 


At this stage, one must note that this integral is divergent. Indeed, ^2 ^2 

for A^O, one has 1/ VI -cos(A) ~ V2/A. Hence the expression Amax = —^— ■ (C.5) 
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One can then immediately perform the integration on from The integration on J'^ is straightforward because of the in¬ 
equation (C.3), so that the diffusion coefficients are given by function (local resonance), and we are left with 




.. ^IT 


dF . 3 ^ fv' 


where the function (J^) is defined as 




1 QE (Ak^)^ 

1 

(mrili)'kl^,. 

1 ” ■^max 




X 


dJ. 


1 


[/m-Q] 


-/d.; 


Jr, J^, J'r, 

X m 


F(J^, J'r, Jr, t)-F(J^, Jr, O^iJ^, J'r, 0 


(D.2) 


After some algebra, the drift coefficients from equation (C.2) are 
given by 




27,! , 

— kn 


(C.7) 


where the function (jjj) is defined as 


f r X d, d 

r, / QE \] 


In -r 

1 ‘57^ 

\rTKcrp\ 


+ m. 


K d 


crj dJ^ 

K 


^+\m\\-Xn 


-^|m;|+lUr max] 


J- m\ Dfmax] 


(C.8) 


where the susceptibility coefficients are generally given by equa¬ 
tion (77). 

The kinetic equation (D.2) is an integro-differential equation 
that governs the evolution of the system as a whole. It describes 
the effects of encounters between any test particle characterized 
by the angle-action coordinates (7^, Jr) and the field particles 
characterized by the (running) angle-action coordinates (7^, 7'). 
Actually, there is no distinction between test and field particles, 
so that they are characterized by the same distribution function 
F{-, t) that evolves in a self-consistent manner, hence the integro- 
differential character of the kinetic equation. This is a character¬ 
istic of the Balescu-Lenard equation describing the evolution of 
the system as a whole. 


These explicit expressions of the diffusion and drift coefficients 
obtained in equations (C.6) and (C.7) allow to estimate in a sim¬ 
ple way the secular flux in the entire action space J-{J^,Jr), 
once we assume that the distribution function is a Schwarzschild 
DF given by equation (25) and that the susceptibility coefficients 
can be approximated by equation (78). 


Appendix D: Relation to other kinetic equations 

The kinetic equation governing the collisional evolution of a sys¬ 
tem of N stars at the order 1/A is the inhomogeneous Balescu- 
Lenard equation (2). This equation conserves the total number of 
stars and the energy, and monotonically increases the Boltzmann 
entropy (H-Theorem). We note that the collisional evolution of 
the system is due to a condition of resonance encapsulated in the 
term In general, this condition can allow for 

local and non local resonances. In the case of tepid discs con¬ 
sidered in the present paper, assuming that only tightly wound 
spirals are sustained by the disc, we justified in equation (72) the 
fact that the resonances are purely local, so that m\-m 2 and 
J^-J^. Furthermore, because of the epicyclic approximation, 
the intrinsic frequencies of the system, given by equation (22), 
depend only on 7^, so that £l-£l(J^). Under these conditions, 
the Balescu-Lenard equation giving the collisional evolution of 
F - F(J^, Jr, t) may be rewritten as 


dF . 2 d [ Y-i 


1 


ttlr 




-fdj'dj; 


X 


dD(j^-j;) 


Jr, j;,J',m-n)l2 

X m 


dF dF 

F{J'^, J'r, Jr, t)-FiJ^, Jr, OjjiJ'^, J'r, 0 


dj' 


. (D.l) 


Appendix D. 1: Fokker-Planck limit 

We can also use this formalism to directly obtain the Fokker- 
Planck equation governing the relaxation of a test star in a bath 
of field stars, assumed to be in a steady state with a distribu¬ 
tion function Fq{J'^, 7'). Proceeding as in Chavanis (2012a), we 
just have to replace in equation (D.2) the distribution function of 
the field particles 7'(7^, 7', t) by the static distribution Fo(J^, J'l.), 
while the time evolving distribution function of the test particle 
is rewritten as ^(7^, Jr, f) for clarity. This heuristic procedure is 
justified in Chavanis (2012a) by an explicit calculation of the 
diffusion and drift coefficients of the Fokker-Planck equation. It 
transforms the integro-differential equation (D.2) into a differen¬ 
tial equation 


dt 



1 

\lDm.m(J,p, Jr, J 4 ,, J'r, m-Sl)\^ 


X in- 


FoiJ^, J'r)^iJ^, Jr, t)-P(J^, Jr, 0^(J^, J'r) 


(D.3) 


which can be interpreted as a Fokker-Planck equation. If we as¬ 
sume that the field particles are at statistical equilibrium (thermal 
bath), described by the Boltzmann distribution 

FoiJ) = ^ (D 4) 


then, using the relation dFojdJ - -(iF^Fl (see the definition of il 
in equation ( 1 )), we can reduce the Fokker-Planck equation (D.3) 
to the form 





+ j0Q(7^)P 


(D.5) 
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with the diffusion coefficient 


DM) = 


4-n^ 




.f4j; 


Fo(J^,J') 


Jr, J'r, m-Sl)Y 


(D.6) 


We note that the friction term in the Fokker-Planck equa¬ 
tion (D.5) is proportional and opposite to the intrinsic frequency 
Jl, and that the friction coefficient ^ satisfies a (generalized) Ein¬ 
stein relation ^ = D/3 (for each resonance) (Chavanis 2012a). 

This Fokker-Planck formalism may have other applications. 
For example, if we consider a system with two species of par¬ 
ticles (e.g. characterized by different masses niA and ms), and 
if species B has reached an equilibrium state with a distribution 
Fo(J'^, J'r), we can use equation (D.3) to describe the relaxation 
of particles of species A due to the encounters with particles 
of species B (but neglecting the encounters between particles 
of species A). This approach is further discussed in Appendix 
F of Chavanis (2013b). On the other hand, for a single species 
system, we may describe the early dynamics of the system as a 
whole with a good approximation by replacing F{J^, J'r, t) in the 
Balescu-Fenard equation (D.2) by the initial distribution func¬ 
tion Fo{J^, J'r), leading again to an equation of the form (D.3). 


Appendix D.2: Other kinetic equations 

We now compare the previous results to other related kinetic 
equations. For spatially homogeneous systems with long-range 
interactions, the Balescu-Fenard equation reads (see Chavanis 
(2012c)): 


dt 


— n{2nym 


dv 


Jdkdv' i 


u(k)^ 


X k- 


\e(k, k-v)]^ 
dF 


6D[k-(v-v')] 


dF 

F{v', t)Mv, t)-F{v, t)—{v', t) 
dv dv 


(D.7) 


where 

e(k,Lo) - l-(2n)‘^u(k) fdv (D.8) 

J kv-u) 

is the dielectric function. For one dimensional systems, it re¬ 
duces to the trivial form 


dt 


— 2n~m 




u(k)^ 


\e(k, kv)[ 


■ 6d(v-v') 


dF dF 

F(v',t)—iv,t)-F(v,t)—iv',t) 
dv dv 


= 0 . 


(D.9) 


Ind > 1, there are always resonances between particles with dif¬ 
ferent velocities, implying that the Balescu-Fenard equation re¬ 
laxes towards the Boltzmann distribution on a timescale of the or¬ 
der Nto, where to is the dynamical time. By contrast, in = 1, the 
resonances become local (in velocity space) and since the term 
in brackets is anti-symmetric with respect to the interchange of 
V and v', the Balescu-Fenard diffusion current vanishes exactly. 
This implies that the relaxation time is larger than Nto, presum¬ 
ably of the order N^to, corresponding to the next order term in 
the expansion of the dynamics in powers of 1 /N. The Fokker- 
Planck equation describing the evolution of a test particle in a 
bath of field particles is discussed in Chavanis (2012c, 2013a). 

The kinetic equation governing the collisional evolution of a 
system of N point vortices in two dimensional hydrodynamics at 


the order 1 /N can be written, in the axisymmetric case, as (see 
Chavanis (2012b)): 


dco 

dt 


1 d r^°° 

= 27r^7- — I dr' r'x(r, r', Q(r, f)) (JD[D(r, t)-Q(r', f)] 
r dr Jo 


^ , 1 (9w 1 doj 

M ,t)-^(r,t)-co(r,f)-^(r ,t) 
r or r or 


(D.IO) 


where w(r, f) is the profile of vorticity, Q(r, t) the profile of angu¬ 
lar velocity, and;^f(r, r', Q.(r, t)) is related to the dressed potential 
of interaction between the point vortices (see Chavanis (2012b) 
for more details). This equation conserves the total number of 
point vortices and the energy, and monotonically increases the 
Boltzmann entropy (H-Theorem). If the profile of angular veloc¬ 
ity is monotonic, the kinetic equation reduces to the form 


dco 

dt 


—2k 


1 (3 0 + 00 1 

'■J--K I dr' r'x{r, r, Q(r, t)) Mr-r') 

r dr Jo |Q'(r, t)\ 


1 dco 1 dai 

cti(r,0-^(r,t)-oj(r,t)-—(r,t) 

y. yjy j./ ^ y. 


r' dr 


= 0 . 


(D.ll) 


For non-monotonic profile of angular velocity, one can have non¬ 
local resonances (i.e. distant collisions between point vortices), 
as studied in Chavanis & Femou (2007). This produces a diffu¬ 
sion current. If the profile of angular velocity is, or becomes, 
monotonic, the resonances are purely local and, since the term 
in brackets is anti-symmetric with respect to the interchange of 
r and r', the diffusion current also vanishes. This implies that 
the relaxation time is larger than Nto as discussed above. The 
Fokker-Planck equation describing the evolution of a test vortex 
in a sea of field vortices is discussed in Chavanis (2012b). 

If we focus on purely local resonances, we note that the in¬ 
homogeneous Balescu-Fenard equation (D.l) is different from 
equations (D.9) and (D.l 1) because it is two-dimensional in Jr 
and J^, and the resonances act only on 7^. Therefore, purely local 
resonances do not yield a zero flux, contrary to equations (D.9) 
and (D.l 1). This really is an effect of the two-dimensionality of 
the system. Indeed, for local resonances, the ID inhomogeneous 
equation also yields a zero flux: 


JL. far 

dt dj\^ Q'(J) J 

\ tr] \ X 


dF dF 

F{J',t)—{J,t)-F{J,t)—{J',t) 


6M-J') 

|2)„,,„(7,7,mD)|2 

= 0 . 


(D.12) 


Appendix E: The Schwarzschild conspiracy 

The Schwarzschild distribution function introduced in equa¬ 
tion (25) and considered in S12 simulation has the specificity 
to be exponential in the /^-direction, so that F is Boltzmannian 
with respect to the Jr variable. It is known {e.g. Chavanis 2012a) 
that the (complete) Boltzmann distribution is the steady state of 
the Balescu-Fenard equation (2). Therefore, we can expect that 
the (partial) exponential behavior of the Schwarzschild distribu¬ 
tion will induce simplifications that we now detail. In analogy 
with equation (10), the flux associated to a given resonance m is 
defined as 




AM)F{J) + DM)m~ 

dj 


= niTm 


(E.l) 


where the non-bold is a scalar. Using shortened notations and 
forgetting numerical prefactors, we may rewrite the drift and dif¬ 
fusion coefficients from equations (79) and (80) under the form 
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where we used the shortened notation 
1 1 


im-m |4-[m-£l] 


(E.2) 


(E.3) 


dJt 


for the term appearing as a prefactor in the expressions (79) 
and (80) of the drift and diffusion coefhcients. The flux can then 
be decomposed as with 


■' in 


(m£l) 








(E.4) 


and 


rrr^ _ 

J m~ 




(mSl) 




\D„ 


T dF , I , 

F{J^r) —(/,") 

OJd) OJd) 


■ (E.5) 


We are interested in the value of the flux at the initial time, 
where the distribution function is given by the Schwarzschild 
distribution. Because of the exponential dependence in 
Jr of the Schwarzschild distribution function, one has 
dFfdJr--(Kl(T^) F. As a result, the radial component (E.4) 
of the flux cancels out and the flux is simply given by equa¬ 
tion (E.5). This coincidence could be called the Schwarzschild 
conspiracy and has important consequences on the properties 
of the collisional diffusion. Indeed, for a tepid disc, one has 
\dF/dJr\ » \dF/dJ^\. Hence one would expect the gradients with 
respect to J, to be the major contributors to the diffusion. When 
considering independently the drift and diffusion coefficients, 
the gradients in dFjdJr dominate the diffusion current. Thus for 
mr + Q, the diffusion-only flux can be approximated by 
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(E. 6 ) 


As the ILR and the OLR have a non-zero nir compared to 
the COR, these resonances should dominate independently the 
drift and diffusion components. However, when considering the 
full flux made of the contributions from the drift and diffu¬ 
sion coefficients, because of the Schwarzschild conspiracy, there 
is a simplification of tfie dominant terms in dFjdJr, so tfiat 
one recovers as in equation (E.5) tfiat only tfie smaller gra¬ 
dients dFjdJ^ remain present. The Schwarzschild conspiracy 
between drift and diffusion will therefore tend to slightly re¬ 
duce the magnitude of the full diffusion flux, so as to moder¬ 
ately slow down the collisional relaxation. More importantly, 
the Schwarzschild conspiracy will favor the COR resonance (ra¬ 
dial migration) over the ILR resonance (y^-heating). One should 
note that the DE which is effectively sampled in S12 is of the 
form E = E(£, y 0 )ocy^exp[-£'/cr^], with q'= Eq/ cr^-1 (Toomre 
1977; Binney & Tremaine 2008). It is only within the epicyclic 
approximation that this DE takes the form of the Schwarschild 
DE from equation (25). As a consequence, in S12 simulation, the 
Schwarzschild conspiracy is not exactly satisfied as observed in 
equation (E.4), but the residual difference driving the secular dif¬ 
fusion is likely to be subdominant, as illustrated in figure E. 1 . 



u 

Eig. E. 1: Contours in action-space (J^, Jr) of the difference between 
the sampled anisotropic DF of S12 simulation F^jOcJ'^ exp[-£/crJ] and 
its Schwarzschild epicyclic approximation Fsch from equation (25). The 
plotted quantity is |F’am“Cschl/F’^“, and contours labels are expressed 
in percentages. 


Appendix F: Temporal frequency selection 

An important feature of the diffusion equation (7) is that the dif¬ 
fusion takes place along specific resonance directions associated 
to the vectors m as discussed in equation (16). Hence being able 
to determine the dominant resonance is crucial in order to esti¬ 
mate the direction of the secular diffusion in action space. The 
temporal frequency associated to a given resonance m in a loca¬ 
tion J of action-space is given hy a> = m £l. Thanks to the expres¬ 
sion (83), one immediately notes that for a Mestel disc, one has 

0 < WiLR < WCOR < WOLR ■ (El) 

In Fouvry & Pichon (2015); Fouvry et al. (2015), we studied the 
same S12 simulation using the WKB limit of the secular diffu¬ 
sion equation, which intends to describe the secular forcing of 
a collisionless self-gravitating system perturbed by an external 
source. An essential assumption of this approach was to consider 
the external perturbation as originating from numerical Poisson 
shot noise, and therefore assume it to be proportional to the local 
active surface density. The autocorrelation of the external pertur¬ 
bation was taken to be equal to 

^|^“f^(w,k„^)(xE,(y^). (F.2) 

One should note that this crude assumption on the noise proper¬ 
ties has no o) dependence, so that all resonances are equally fa¬ 
vored by the Poisson shot noise and are perturbed similarly what¬ 
ever their associated intrinsic frequencies a>-m £l. This ad hoc 
and simple noise approximation is one of the limitations of the 
formalism presented in Fouvry & Pichon (2015); Fouvry etal. 
(2015). In contrast, in the WKB Balescu-Lenard equation de¬ 
scribed in this paper, this preferential selection of the resonances 
based on their intrinsic frequency is naturally present. Indeed, 
one can note in the expression (E.5) of the flux associated to 
a resonance m, the presence of the prefactor llim-il)' which 
arose in equation (67) when handling the resonance condition. 
For the Mestel disc, whose intrinsic frequencies are given by 
equation (83), this term can be straightforwardly computed and 
reads 


(m-Q)' \m^+ sj2mr\ \d0.ri,ldJ,p\ ' 
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Comparing the ILR resonance to the OLR and COR, one imme¬ 
diately obtains that 

(wtQLR-^^)' _ 2-hV2 ^ 2 g (wtcoR-^)' _ 2 .^3 4 

(miLR-n)' 2 -V 2 ’ ’ (»iiLR'Si)' 2 -V 2 

(F.4) 

Hence because the ILR resonance is associated to lower intrin¬ 
sic temporal frequency at-m il, the resonant factor l/(»i Q)' 
naturally tends to favor the ILR resonance with respect to the 
OLR and COR, and therefore performs natively a temporal 
frequency biasing which was absent from the ad hoc assump¬ 
tion of equation (F.2) describing the external forcing considered 
in Fouvry & Pichon (2015); Fouvry et al. (2015). 

One can even be more specific when comparing the ILR 
and OLR resonances. The only difference between these two 
resonances is the sign of m^. The expression (62) of the am¬ 
plification eigenvalues shows that its value only depends on s^, 
so that /1 ilr = /Iolr- The expression of the susceptibility coeffi¬ 
cients from equation (77) is also independent of the sign of nir 
so that 1 /IDilrI^ = 1 /|£)olrP- Hence, when considering the flux 
Tm given by equation (E.5), one notes that between the ILR and 
OLR resonances, Tm only changes through the factor l/(m fl)'. 

Thanks to equation (F.4), one immediately obtains 

(F.5) 

TOLR 

The secular diffusion flux associated to the OLR resonance is 
therefore always much smaller than the one associated to the ILR 
resonance, because of this effect of temporal frequency biasing. 

As a conclusion, the temporal frequency selection effect de¬ 
scribed in equation (F.4) will tend to favor the ILR resonance 
because it is associated to a smaller intrinsic frequency. How¬ 
ever, one can note from figure 8 that the COR resonance is al¬ 
ways more amplified than the ILR resonance. Finally, the cru¬ 
cial remark is to note that the susceptibility coefficients from 
equation (77) involve Bessel functions J'm,, which are such that 
= 1 if nir-Q, or = 0 otherwise. As a consequence, 
close to the 7^ = 0 axis, the COR resonance will always tend to 
become the dominant resonance. There is therefore a non triv¬ 
ial arbitration between these opposite effects when considering 
the respective contributions of the various resonances to the full 
secular diffusion flux. 
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